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A  .   I  NT  P.  ODU  C  T I  ON 

The   procedure   for   the   design  of   lowpass,   highpass . 

bandpass    and    bandstop    digital  filters    from   lowpass 

prototypes  is   well  established,  and  is  based  en  transferring 


ths"    prototype    frequency    response 


to 


r.  e  e ' 


desired 


C.»=.^  c  • 


This    t r ans f orTTH t ion    can   be    accoir.Dlished    in 


either    the      analog       (s)       dorrain      or      the      digital       fz)    dorr.ain 
(Figure    1 ^  . 


lowr-ass       ' 
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(DIGITAL  TRANSFORMATIONS )=^  digital' 
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Figure  1.   Two  methods  for  digital  filter  design. 

Note  that  each  of  the  horizontal  arrows  actually 
represent  four  different  transformations: 
1^  lov^'oass  prototype  to  lowpass  filter 
2)  lowDass  r>rototype  to  hiohpass  filter 
3^  lovrpass  prototype  to  bandpass  filter 
4^  lowpass  prototype  to  bandstop  filter 


Therefore.  the  two  classical  design  techniques.  analog- 
diaital  and  digital-digital  design,  are  based  on  doing  the 
frequency  response  transformations  in  their  respective  analoa 
or  digital  domains.  Since  the  theory  for  the  standard  filter 
types  (Butterworth ,  Chebyshev,  and  elliptic)  is  developed  in 
the  analog  or  's'  domain,  the  two  methods  for  designing 
digital  filters  differ  basically  in  the  sequence  of  the 
design  steps.  In  the  analog-digital  method  the  analog 
prototype  is  transformed  and  then  digitized,  while  in  the 
digital-digital  method  the  analog  prototype  is  digitized  and 
then  transformed. 

Each  of  the  transformations  (arrows  in  Figure  1)  involves 
substituting  a  function  of  's'  or  'z'  for  the  appropriate 
variable  in  the  filter  transfer  function.  The  purpose  of 
this  thesis  is  to  develop  efficient  algorithms  to  carry  out 
these  transformations  and  then  to  implement  these  algorithms 
in  a  computer  program  to  aid  in  the  design  of  digital 
filters . 

Com.puter-aided  desian  of  digital  filters  offers  a 
significant  advantage  since  it  allows  the  digital  filter 
coefficients  to  be  calculated  quickly  and  accurately.  The 
accuracy  of  the  filter  coefficients  is  an  especially 
significant  point  since  the  frequency  response  of  the  final 
filter  is  very  sensitive  to  sm.all  inaccuracies  in  the 
calculated  coefficients.  As  an  example,  compare  the 
freauency  response  for  the  filters  in  Figure  2  and  3.   Fiaure 


2  shows  the  frequency  response  magnitude  for  a  1/2-dB  ripple 
Chebyshev  bandpass   filter,  and   Figure  3  shows  the  frequency 
response  when  one  of  the   coefficients   has   been   changed  by 
1/lOOOth.       .-    . -  .      
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Figure   3. 
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Same  Filter  With  Denominator  Coefficient  of 
z  changed  to  -.67594 


B.   GENERAL  DISCUSSION 

Each  of  the  transformations  involves  substitutina  a  ratio 
of  two  polynomials  of  s  or  z  into  the  filter  transfer 
function  which  is  also  a  ratio  of  two  polynomials  in   s  or  z , 


namely 


H(y)  = 


Y(y) 
U(y) 


p(x) 


y= 


(1.1) 


a(x) 


The  variables  'x'  and  'y'  can  be  either  's'  or  'z'  depending 
on  the  transf orm.ation  to  be  performed.  Therefore,  we  have 
either 


an  +  ai  y  +  a2y2  +  ...  +  arry" 


H(y) 


bn  +  bi  y  -t-  bz  y2  + 


+  bn  y" 


p{x) 


(1.2) 


or 


a  (x! 


p(x)      (p(x))2 
ao  +  ai  +  32 + 


a(x) 


(a(x) )2 


H(x)  = 


p(x)      (p(x))2 
bo  +  bi  +  bi + 


a(x) 


(a(x) )2 


am 


(D(X) ) 


(a{x) ) 


+bn 


(p(x) )" 


(a(x)  )'■ 


1.3) 


In  general,  the  order  of  the  numerator  'm' ,  is  less  than  the 
order  of  the  denominator  'n'  (except  for  the  transformations 
in   the   diaital   domain   where   the   result   of  the  bilinear 


transf ormiation  is   to   make   n=m 


We   can   generalize  the 


results,  however,  by  assuming  that  the  order  or  the  numerator 
polynomial  YCx)  is  equal  to  the  order  of  the  denominator 
polynomial   U(x)    and   making   the   additional   coefficients 

am* 1, am* 2 an  equal   to  zero.    Then   multiplying  both  the 

numerator  and  denom.inator  of  the  above  equation  by  (g(x))" 
gives 


ao(a(x))"  +  ai  p  (x)  (g  (x)  ) "  -  1  +  .  .  .  +  an(p(x))'' 

H(x)=  (1.4) 

bo(a(x))"  -t-  b]  p(x)  (g(x)  )f'- i  +  .  .  .  +  bn(p(x))° 


or.  in  a  different  form. 


i=n 

I    ai  p(x) 1 g (x)°- ' 
i  =  0 

H(x)=  .  (1.5) 

i=n 

Z   b,  p  (x)  ^  g  (x)  '^  -  ^ 
i  =  0 


Now  with  the  numerator  and  denominator  of  the   same  formi, 

we  can  use  the  same  transformation  for  both  the  numerator  and 

denom.inator   polynom.ials   of   the   filter   transfer  function, 

namely, 

(1.6) 
i=n 
a'(x)=  an'*  ai'x-t-  a2'x2+.  .  .+  ar'x""  =  I  a,  p(x)*g(x)°-» 

i  =  0 

In  Equation  1.6.  the  ai  are  the  coefficients  of  the  old 
numerator  or  denominator  polynomials  and  ai  '  are  the  new, 
transformed  numerator  or  denom.inator  polynomial  coefficients. 
Since  the  new  coefficients  ai  ' .  are   a  linear   combination  of 


the   old   coefficients   at  ,   and  the  coefficients  of  p(x)  and 
a(x),  this  suggests  the  use  of  a  matrix  transformation; 


where 


a"  =  aA(n) 


a     is  the  coefficient  row  vector  for  the  old 
numerator  or  denor.inator  polynomial  of 
lenath=n+l . 


(1.7) 


a  =  r  an  a  i  a? 


.  an  1 


a'    is  the  coefficient  row  vector  for  the  new 
numerator  or  denominator  polynomial  of 
lenath=n ' +1 . 


a'  =  [a'o  a'l  a'? 


an  •] 


A(n)  is  the  transf orm.ation  matrix  with  dimensions 
(n+l ) X (n ' +1 ) .  The  elements  of  A{n)  depend  upon  the 
coefficients  of  p(x)  and  a(x).   Note  that  if  the 
rows  and  columns  of  A(n)  are  num.bered  0  to  n'  ,  then 


i  =  n 
'  =   Z    ai  A  (n )  ,  1 
i  =  0 


1.8) 


(where  A(n)i  i  denotes  the  element  in  row  i, 
column  j  of  m.atrix  A(n).) 

n     IS  the  order  of  the  transfer  function  before  the 

transformiation 

n'    is  the  order  of  the  transfer  function  after  the 
transf  orm.ation .   Note  that  n'  equals  n  times  the 
highest  order  of  p(x)  and  g(x). 

For  example,  when  n=l  and   n'=2,   p(x)   and   g(x)   are  second 

order  and  Eauation  1.7  becomes, 


la. 


a-  I   IA«'2)oo    Af2)o,    A(2)n2 


A  (  2  )  :  f,    A  (  2  )  !  1    A  (  2  )  1  !• 


=  laO'   al'   a2' 


(1.9) 


In  the  case  of  the  lowpass  digital,  highpass  digital,  and 
bilinear  transformations,  the  two  functions  p{x)  and  g(x)  are 
first  order.  Professor  P.H.  Moose  of  the  Naval  Postgraduate 
School  of  Monterey  has  developed  the  rriatrix  A(n)  for  these 
two  cases  by  expanding  Equation  1.6  using  the  binon^ial 
theorem  (Reference  1).  This  method  cannot  be  used  for 
transformations  involving  higher  order  polynomials  such  as 
the  bandpass  and  bandstop  transformations  and  it  is  difficult 
to  develop  because  intricate  manipulation  of  the  binom.ial 
coefficients  are  involved.  The  binomial  expansion  approach, 
therefore,  was  not  pursued. 

A  general  method  has  been  found  that  develops  the 
required  transformation  matrices  in  an  iterative  manner.  In 
other  words ,  the  matrix  used  to  transform  a  numerator  or 
denominator  filter  polynomial  of  order  n  is  developed  from 
the  matrix  for  the  polynomial  of  order  n-1.  This  approach  is 
easy  to  im.plement  in  a  computer  program,  and  can  be  used  for 
any  of  the  polynomial  substitutions  needed.  The  next  chapter 
is  concerned  with  developing  the  transformation  matrices  for 
each  specific  transformation.  However,  as  an  introduction  to 
tne  method  used,  consider  the  generic  case  where  both  pix) 
and  g(x)  are  second-order  polynomials.  Since  all  of  the 
analog  and  digital  transformations  involve  substitutions  with 
polynomials  of  second  order  or  less,  this  second-order 
aeneric  case  will  be  all  that  is  needed  to  develop  the 
specific  transf orm.ations  later. 


With  p(x)=ax2  +  bx  +  c  ,  g(x)  =  dx^ +  ex  +  f  .  and 
substituting  into   Equation  1.6   for  the   first  order  case  of 
n=l  (remembering  that  with  p{x)  and  g(x)  of  second  order,  the 
resulting  polynomiial  will  be  second  order),  we  get 

a'  (x)  =  ao  '  +  ai  'X  +  az  '  x^  =  an  (dx^  +  ex  +  f )  +  ai  (ax^+bx+c) 

(1.10) 

Using  Equation  1.8  or  1.9  it  is  easy  to  pick  out  the  matrix 
coefficients  A(l)ii  (i=0,l;  j=0,l,2)  for  the  first  order 
Generic  case.   The  coefficients  are, 


Ad)  = 


f 

c 


e 
b 


d 

a 


(1.11) 


Now  with  n=2  (n'=4)  substituting  into  Equation  1.6  gives. 


'  (x)  =  ao  '  +  ai  'x  +  a2  'x'  +  as  '  x^'  +  a^  '  x' 


(1.12) 


=  ao  (dx' +ex+f ) ^  +  ai  ( ax^ +bx+c ) (dx^ +ex+f )  +a2  ( ax^ +bx+c ) ^ 


and  the  matrix  A(2)  is; 


If2 

2ef 

2df+e2 

2de 

d2  1 

A(2)  = 

jfc 

ec+bf 

cd+f a+eb 

db+ea 

da  i 

2bc 

2ac-»-b' 

2ab 

a  i 
1 

(1.13) 


To  see  how  this  matrix  is  derived  from  the  previous  one, 
consider  the  first  row  of  the  matrices  in  Equations  1.11  and 
1.13.  For  both  matrices,  this  row  is  multiplied  by  the  ao 
coefficient  in  the  old   polynomial.    For  n=l   we  have  aop(x) 

8 


and  we  pick  out  the  coefficients  of  the  powers  of  x  (x-  )  to 
get  the  matrix  elements  A{n)oi.  For  n=2 ,  we  have  ao  (p(x))^ 
and  we  do  the  same  thina  only  now  we  have  the  additional 
factor  of  p(x).  The  matrix  elements  can  be  calculated  by 
hand  in  this  fashion  but  it  requires  a  lot  of  algebra.  A 
simpler  method  is  to  use  the  matrix  elements  of  the  previous 
matrix  and  get  the  appropriate  elements  by  convolving 
coefficients.  For  example,  consider  the  sequence  of  elements 
in  row  0  (the  top  row)  of  Ad)  ,  namely 


(1.14) 


To  get   the  elements   for  row   0  of  A(2) ,  we  need  to  multiply 
(convolve  coefficients)  by   another   p(x).     Therefore,  with 
rows  numbered  C  to  n,  and  columns  number  0  to  n' ,  we  have  the 
fol lowing. 

For  row  0  colum.n  0,  the  appropriate  matrix  element  is 


A(2)oo  =  d    e    f  =  f2 

fed  (1.15 


For  row  0  column  1  ,  the  appropriate  matrix  element  is 


A (2)o .  =      d    e    f  =  2fe 

fed  (1.16) 


For  row  0  coluT.n  2,  the  appropriate  matrix  element  is 


A(2W,  2  =  d        d    e    f  =  2fd+e2 

fed  (1.17) 


For  row  0  column  3,  the  appropriate  matrix  element  is 

A{2)o3=  d    e    f       =  2de 

fed  •  (1.18) 


For  row  0  column  4 .  the  appropriate  matrix  element  is 


A  ( 2 )  0  4  =  d    e    f   =  d2 

fed  (1.19) 


This  convolution  of  the  coefficients  of  the  old  matrix  with 
the  coefficients  of  p(x)  can  be  carried  out  for  every  row  of 
the  old  matrix  to  obtain  all  the  rows  of  the  new  matrix 
except  the  last  row.  In  the  case  of  the  last  row,  we  are 
multiplying  by  an  additional  g(x)  rather  than  an  additional 
p(x)  and  so  must  obtain  the  last  row  of  the  new  matrix  by 
convolving  the  coefficients  of  the  last  row  of  the  old  matrix 
with  the  coefficients  of  g(x). 
In  sur.mary  we  have , 

(i)  for  row  0  to  row  n-1: 


•row  i  of  A(n)  =  row  i  of  A(n-l)  *  coefficients  of  p(x) 

(1.20) 


( ii )  for  row  n 


-row  n  of  A(n)  =  row  n-1  of  A(n-l)  *  coefficients  of  g(x) 

(1.21) 


where  *  indicates  linear  convolution  of  the  coefficients. 

To  see  how  this  method  of  convolving  coefficients  works. 
consider  multiplying  two  polynomials: 

fl(x)=f2(x)f3(x)  (1.22^ 

10 


Since  'x'  is  the  independent  variable,  no  matter  what  letter 
we  use,  the  polynomials  will  remain  the  same.  Substitut inc 
' z" ^  *  for  'x'  in  the  above  equation  gives  us. 


fl  (2- M  =  f2(z-i  )f3  (2-^  ) 


(1.23) 


and  we  now  have  the  familiar  form  where  we  can  either 
multiply  the  two  functions  of  '2'  or  convolve  the  sequences 
found  by  taking  the  inverse  2-transform.  For  f2(2-M  and 
f 3 { 2- M  we  find 


f2{2-i)  =  ao  +  ai2-i  +322-2  +  . 


f  3  (2-  M  =  bo  +  b]  2-  ^  +  b2  2-2  + 


+  an  2-  0 


+  b»2-«' 


(1.24) 


It  is  now  clear  that  we  can  find  the  coefficients  of  the 
product  fl(x)  by  convolving  the  sequence  Ian  i  with  the 
sequence  I bm I  to  produce  the  sequence  1 Ca  I  where. 


an  I  =  I  ao  ,  ai  ,  az  ,  33  , 


I  bm  I  =  i  bo  ,  bi  ,  b2  ,  b^  , 


and  with  a  •=  m  +  n 


i  Cq  I   =   f  Co  ,  C1  ,  C2  ,  C-?  , 


,  an  I 

,  bir  I 

.  Cq  > 


(1.25) 


Then   by   taking   the   2-transform   of   the  sequence  I Co  }  and 
substitutina  'x'  for  '2-' '  we  aet , 


fl(x)   =  C(-   +  C)X  + 


+   Cq  X*' 


(1.26) 


11 


We  now  have  a  simple  straight-forward  method  of 
developing  the  necessary  transformation  matrices  for  digital 
filter  design  by  picking  out  the  first-order  transformation 
matrix  using  Equations  1.6  and  1.11,  and  then  using  the 
convolution  of  coefficients  method  iteratively  on  the  rows  of 
each  successive  transformation  matrix  with  Equations  1.20  and 
1.21  until  the  appropriate  matrix  A(n)  has  been  generated. 
The  next  chapter  will  develop  the  transformation  matrices  for 
each  of  the  following  transformations: 
bilinear  transformation 

digital  lowpass  prototype  to  digital  lowpass  filter 
digital  lowpass  prototype  to  digital  highpass  filter 
digital  lowpass  prototype  to  digital  bandpass  filter 
digital  lowpass  prototype  to  digital  bandstop  filter 
-   analoo  lowpass  prototype  to  analog  lowpass  filter 
analog  lowpass  prototype  to  analog  highpass  filter 
analog  lowpass  prototype  to  analog  bandpass  filter 
analog  lowpass  prototype  to  analog  bandstop  filter. 
We  will  see  that   sone  of   these  transformations   are  sim.ilar 
and    consequently,     we    will     not    have    to    develop 
seoarate  transformation  matrices  for  each  case. 
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II.   pEyELOPMENT  OF  SPECIFIC  TRANSFORMATIONS 

A.  INTRODUCTION 

The  individual  transformations  are  discussed  in  the 
following  sections.  The  specific  form  of  the  transformations 
used  in  each  case  are  from  Reference  2.  Since  the  final 
result  will  be  to  develop  a  computer  program  to  assist  in  the 
calculation  of  digital  filter  coefficients,  confining  the 
specific  transformations  to  the  form  in  Reference  2,  will 
allow  the  computer  program  to  be  used  in  conjunction  with 
Reference  2  as  an  aid  in  designing  digital  filters. 

B.  BILINEAR  TRANSFORMATION 

The  bilinear  transformation  may  be  used  to  transform  an 
analog   transfer    function,   H(s) ,   to   a   digital   transfer 

function  K(z).    The   transformation   from   H{s)   to   H(z)  is 
described  by 


H{z^  =  H(z)  I 

i     z-1  (2.1) 

!s=  

z  +  1 


Referring  to  Equation  1.1,  p(x)=z-l  and  g(x)=z+l.  Since  both 
of  the  transf orm^ation  polynomials  p(x)  and  g(x)  are  first 
order,  the  bilinear  transformation  will  not  change  the  order 
of  the  transfer  function.  Therefore,  for  an  analog  transfer 
function  of  order  'n' ,  the  transformation  matrix  A(n)  will  be 
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(n+l)  by  (n+1).  From  Equation  1.6,  the  numerator  or 
denominator  polynomial  of  the  transfer  function  H(z)  is  given 
by 


a  '  (z)=  ao  '  +  ai  •  z  + 


i=n 
+  an  "z"  =   I   ai  (z-l)i  (z+1)"-* 
i=0 


(2.2) 


where  the  unprimed  coefficients,  as  ,  come  from  the  analog 
transfer  function's  numerator  or  denominator  polynomial.  In 
matrix  form,  with  n=l  and  referring  to  Equation  1.11,  the 
first  order  matrix  A(l)  for  the  bilinear  transformation  is 


A(l)  = 


1 
■1 


1 

1 


(2.3) 


Usina  the  convolution  of  coefficients  method  developed  in 
Chapter  1.  the  transf ormiation  matrix  A(n+1)  can  be  developed 
from  the  matrix  A(n).  The  results  for  the  bilinear 
transformation  matrices  Ad)  through  A(5)  are  given  in  Table 
1. 
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TABLE  1 
FIRST-ORDER  TO  FIFTH-ORDER  BILINEAR  TRANSFORMATION  MATRICES 


Ad)     =        I     1 
i 

!-l 
L_ 

A{2)    =       I    1 
-1 


1 
1 

2 

0 

•2 


A(3)     = 


A{4)     = 


1 
-1 

1 
-1 

1 

-1 
1 

-1 
1 


3 
■1 
•1 

3 

4 

•2 

0 

2 

■4 


3 

1 
■1 
■3 

6 
0 
•2 
0 
6 


1 
1 

1 
1 

4 

2 

0 

-2 

-4 
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i 

1  i 
i 

i 
li 

i 

1  i 


A(5)     = 


1 
•1 

1 
■1 

1 
•1 


10 


3 
1 
1 
3 
5 


2 

■2 
2 
2 

■10 


10 
2 

-2 

-2 
2 

10 


5 

3 

1 

■1 

■3 


1 
1 

1 
1 
1 
1 
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Referring  to   Table  1  and  keeping  in  mind  the  convolution 

of   coefficients   method  that   produced  the   matrices,   the 

following   relations   are  apparent   and  will   be   useful  in 
computer  implementation. 

A(n)oi  (the  first  row)  =  (°)     for  j=0.1 n        (2.4) 

where  (  ?  )  indicates  the  binomial  coefficient 

A(n)io  (the  first  column)  =  (-1)^   for  i  =  0 , 1 n      (2.5) 

A(n)in  (the  last  column)  =  1     for  i=0,l n         (2.6) 

A(n)(n-i)i  =  (-1) '' -  J  A  (n)  1  )  for  i  =  0  , 1 ...,  integer  (n/2  ) -H  (2.7) 


A(n).i  =  A(n-l),i  +  A{n-l)i(i-i>   for  i=0 . 1 , . . . , n-1     (2.8) 

j=l , 2 , . . . , n-1 


Notice  that  the  values  of  '  the  matrix  elements  for  the 
bilinear  transformation  are  whole  numbers.  This  is  important 
from,  a  numerical  accuracy  standpoint  since  whole  numbers  can 
be  expressed  precisely  in  floating  point  form.  This  means 
that  the  bilinear  transformation  will  cause  a  loss  of 
accuracy  of  the  original  coefficients  only  because  of  the  n+1 
multiplications  and  additions  used  to  calculate  each  new 
digital  coefficient. 

C.   LOWPASS  TO  LOWPASS  DIGITAL  TRANSFORMATION 

The  lowpass  to  lowpass  digital  transformation  is  used  to 
transform  a  lowpass  prototype  digital  filter  H(z)  into  a 
lowpass  digital  filter   H{z).     The   polynomial  substitution 
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used  is 

H(2)=  Hp  (z) I 

z-a 

|z=  (2.9) 

1-az 

where  the   subscript  p  indicates  the  prototype  digital  filter 

transfer  function.     The   design   constant   a   is  determined 

[References  2  and  3]  by 

sin(er  /2  -  Sc  '  /2) 


a  = 


sinOc  /2  -  er  "/2)  (2.10) 

with  Gr  =  the  prototype  critical  frequency  (usually  pi/2)  and 
©r  '=  the  desired  critical  frequency  of  the  lowpass  filter. 
Again,  referring  to  Equation  1.1  with  p(x)=  z-a  and  g(z)=  1- 
az ,  the  lowpass  digital  transformation  matrix  will  not  change 
the  order  of  the  transfer  function  and  the  transformation 
matrix  A(n)  will  be  (n+1)  by  (n+1).  Using  Equations  1.6  and 
1.11  as  before,  the  first  order  transformation  matrix  A(l) 
for  the  lowpass  digital  transformation  is: 


i  1    -al 
A(l)  =  i         I  (2.11) 

i-a     1  I 


Once  again  using  the  convolution  of  coefficients  method,  the 
higher-order  transf orrnation  matrixes  A(n)  are  developed  from 
the  lower-order  matrix  A(n-l) .  The  results  are  shown  in 
Table  2  for  the  first-order  to  fifth-order  transformation 
matrices . 
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TABLE    2 
FIRST-ORDER    TO    FIFTH-ORDER    LOWPASS    DIGITAL 
TRANSFORMATION    MATRICES 


A(l)= I     1  -a 
i 

A(2)=  I     1  -2a           a2  i 

I  ! 

l-a  (l+a2  )       -a  I 

i  I 

ia^  -2a              li 

L_  _l 

A(3)  =  |    1  -3a                3a2               -a^ 

i    -a  (l+2a2 )        (-2a-a3 )         a^ 

i 

Ia2  (-2a-a3  )     (l+2a2 ) 


-a 


I 


i-c3  3c2                  -3a                    1     ! 

i_  _l 

'  "^                           i 

l-a  (l-^3cr2)             (_3c_3c,3  \           (302-^0"^          -a^  i 

1  i 

1^2  (-2~-2a'')        (l-*-4c:2-»-a''^        (-2a-2-^)          a2i 

1  ~                        I 

\-ry3  (3a2-i.a4)           ^-3a-3a3^           ('i+'^a2x             -a     ' 

i 

'a*'  -4o3              6a2                          -4a                         1  ! 

L_  _i 

A(5i='     1  -5a                    10a2                          -lOa^                       So"             -a'  ' 

'  i 

i-a  M +?a2  ^           (_4^_6(,.3  )                (6a2+aM           (-4o3-a=)       a4  I 

1  i 

Ia2  {-2a-3aM     (l  +  6a2 +3aM        {-3a-6a3-a')     (3a2+2a'')    -aM 

i  j 

!-a  (3a2+2aM     (-3a-6a2-a5)     (l  +  6a2+3a'')        (-2a-3a3)       a2  I 

i  I 

la^  (-4a3-a5)     (6a2+4a'')                (-4a-6a3)           (l+'4a2  )          -a     i 


-a^ 


Sa"  -10a3  10a2  -5a 
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Referring   to   Table   2,   the   following    relations   are 
observed  and  will  be  used  in  the  computer  implementation. 


A(n)oj   (the  first  row)  = 


(  "  )  (-a) J   for  j  =  0,l, 


,n   (2.12) 


A(n)io   (the  first  column)  =  (-a)*   for  i-0,l,...,n     (2.13) 


A(n)i)   =  A (n) ( n -1  )  < n- )  )   for  i,j=o,l, 


n 


(2.14) 


A(n)ii   =  A(n-l) (  1  -  1  ) <  1  - )  )  +  (-a) A (n-1 ) <  j - i  >  i 

for  i=l , 2 ,    , n-1 
j  =  l  ,2 n 


(2.15) 


Note  that  the  elements  of  the  lowpass  digital 
transformation  matrices  are  not  whole  numbers  as  in  the 
bilinear  transformation  case,  but  are  polynomials  of  the 
design  constant  a.  Numerical  error  is  introduced  in  this 
case  by  both  the  n+1  additions  or  multiplications,  and  the 
error  in  calculating  the  matrix  elements  that  multiply  the 
original  coefficients. 

D.   LOWPASS  TO  HIGHPASS  DIGITAL  TRANSFORMATION 

The  lowpass  to  highpass  digital  transformation  is  used  to 
transforiTi  a  lowpass  prototype  to  a  highpass  filter.  The 
polynomial  substitution  for  this  transformation  is: 


H(z)  =  Hp  (z 


z-a 


z  =  - 


1-az 


(2.16) 


Since  this  transformation  is  the  same  as  that  for  the  lowpass 
to  lowpass  digital  transformation  except  for   the  minus  sign, 
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we  can  make  the  substitution  of  ' -z'  for  '2*  in  the  original 
prototype  transfer  '  function  and  then  use  the  saFie 
transformation  matrices  as  the  lowpass  to  lowpass  case.  For 
the  lowpass  to  highpass  transformation,  the  design  constant  a 
is  (References  2  and  3) 


cos (0c /2  -  ec  ' /2) 

a  =  (2.17) 

cos  (Sr  /2  +  Gr  '  /2) 


with  Gc  -  the  prototype  critical  frequency  (usually  pi/2)  and 
Gc '  -  the  desired  critical  frequency  of  the  highpass  filter. 

E.   LOWPASS  TO  BANDSTOP  DIGITAL  TRANSFORMATION 

The  polynomial  substitution  which   transforms   a  digital 

lowpass  prototype  to  a  bandstop  digital  filter  is  (References 

2  ar-d  3 )  : 

2a       1-k 

H(z)  =  Hr  (z) I     z^  -  z  +  

I  1+k      1+k  (2.18) 

I2  =  

2a       1-k 

1  -   z  +  z2 

1+k      l+k 


with 


cos  (0,.  '  /2  +  Gi  •  /2) 

a  =  (2.19) 

cos (Gu  ■ /2  -  Gi  ' /2) 


ana 


K  =  tan(Gc /2) tan(Gu  '/2-0i'/2)  (2.20) 

with   once   again,   Gc   =   the   prototype   critical  frequency 
(usually   pi/2),   0n  '   =   the   desired  upper  critical  digital 
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frequency,   and   Gi  '   =   the   desired   low   critical   digital 

frequency. 

By  letting   b=  2a/(k+l)   and  c= ( 1-k ) / ( 1+k ) .  Equation  2.18  can 

be  simplified  to 


H(z)  =  Hp  (z) 


2  = 


z2  -  bz  +  c 


1  -  bz  +  cz2 


(2.21) 


Again,   using   Equations   1.6    and    1.11    the   first-order 
transf orination  matrix  is 


A(l)  =  II    -b 


I  ^ 


-:b 


(2.22) 


Mote   th?t   since   the   subst iti:tior.   polynorrials   are  second 
order,  thi2   transf orratior   will   double   the   order   of  the 


C-^  A  f~,  A  y-  ■s  "^  •-»»-~*-^*-  'r»^  ^ 


^  A    1    *-  c->- 


1'^er^f ore  ths  transf or'^aticr* 
Tratrix  A'n^  will  be  (n+1^  by  (2n+l).  The  higher  order 
transformation  matricies  are  again  developed  recursively  and 
the  results  for  the  first-order  to  third-order 
transformations  are  aiven  in  Table  3. 
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A(l)     = 


A(2)  = 


!c2 
1_ 


TABLE    3 

FIRST-ORDER    TO    THIRD-ORDER    BANDSTOP 

DIGITAL    TRANSFORMATION    MATRICES 


■b  c 

■b  ] 

-2b  (b2+2c)  -2bc         c2  i 

i 

(-b-bc)     {l+b2+c2)     (-b-bc)       c     j 
-2bc  (b2+2c)  -2c  1     I 


_l 


A(3)    =      1  -3c  {3b2+3c)                      (-b^-Gbc) 

Ic  (-2bc-b)  (l+2b2+bc+2c)        ( -b^ -2b-2bc-2bc2 ) 

J 

ic2  (-2bc-bc2 )  {c3 +2b2 c+2c+b2 ) 


-3c2b 


(3b2  c  +  3c2  ) 


(Redundant  entries  not  shown  for  A{3)  due  to  space 
lirr.itations .   Entries  not  shown  can  be  found  using  Equation 
2.25) 


Notice  the  increased  amount  of  calculation  needed  to 
deterir.ine  the  ifiatrix  elements  for  the  bandstop  case.  Each 
matrix  element  is  a  function  of  'b'  and  'c'  which  are 
themselves  functions  of  the  design  constants  a  and  k.  Table 
3  only  gives  the  transformation  matrices  for  first-  to  third- 
order.  However.  to  illustrate  the  amount  of  calculation 
required,  the  matrix  element  A(5)25  is 

A(5)2n  =  (-12bc2 -6bc= -8b3  c-6bc-6b3 -3b-3bc'' -6b3c2-b=  )    (2.23) 


Since   each   coefficient   in   the   final   filter   is  found  by 
multiplying   n+1   matrix   elements    by    the    n+1   original 
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coefficients  and  then  adding  the  n+1  products,  the  accuracy 
of  the  design  constants  and  the  original  coefficients  can  be 
critical  in  the  calculation  of  the  final  filter  coefficients. 
The  symmetry  relations  for  the  bandstop  transformation 
are : 


A(n)io   (the  first  column)  =  c*   for  i=0,l,...,n 


(2.24) 


A  { n )  i  1  = 


A{n)(n-i)(2n-i)   for  i=0,l,...,n 

j  =  0.1 2n 


(2.25) 


A(n)i  1  =  A(n-l) t   i 


-  bA (n-1 ) i  (  1  -  1  )  +  cA(  n-1 ) j  (  1  -  2 ) 

for  i  =  0  ,1 n-1 

j  =  0  ,1 2n 


(2.26) 


F.   LOWPASS  TO  BANDPASS  DIGITAL  TRANSFORMATION 

Again,  the   polynomial  substitution   which  transforms  the 
lowpass  prototype  to  a  bandpass  filter  is 


H(z)  =  Hp  (z 


2ak      k-1 
:2  -  z  +  


k+1 


k+1 


(2.27) 


z  =  - 


1  - 


2ak 


k+1 


z  + 


k-1 


k  +  1 


where  in  this  case  k  =  tan  (9/2 )  cot  (6.,  ' /2  -  Oi  '  /2)  and  a  is 
found  from  the  sam.e  expression  as  the  bandstop  case.  With 
the  simplification  of  b=2ak/(k+l)  and  c= (k-1 ) / (k+1 ) ,  Equation 
2.27  becomes 
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H(z)  =  Hp  (z) 


2  =  - 


z2  -  bz  +  c 


1  -  bz  +  cz2 


(2.28) 


Since  this  is  the  same  form  as  the  bandstop 
transformation,  we  can  use  the  matrices  previously  developed 
if  '-z'  is  first  substituted  for  'z'  as  was  done  in  the 
lowpass  and  highpass  situation. 

G.   LOWPASS  TO  LOWPASS  ANALOG  TRANSFORMATION 

The  lowpass  to  lowpass  analog  transformation  is  used  to 
transform  an  analog  lowpass  prototype  to  a  lowpass  analog 
filter.   The  polynomial  substitution  is  a  simple  one,  namely 


H(s)  =  Hp  (s) 


(2.29 


S  =  S/W( 


where  Wr  is  the  prewarped  critical  frequency  (wr=tanOr/2) 

and  0r  is  the  desired  digital  critical  frequency 

(Gr =2 (pi ) f c /f s  ) .   This  transformation  matrix  is  found  to  be 


1      0 

0       1/Wr 


A  ( n  ;  =  I  C 


lo 


0 
0 


1/Wc 


0 
0 
0 


l/WrM 


or  A(n'  =  diaa(  1,  l/Wc   ,1/Wc 


,l/wc">  throughout. 

(2.30) 


24 


Note  that   in   this   case,  each   new   filter  coefficient 

depends   only   on   multiplying  the   old  coefficient  of  s'  by 

(wi  )"' .   Therefore,  any  loss  of  accuracy  due  to  the 

transf orjTiation  is  due   to   this  one   multiplication   and  the 

error  introduced  in  calculation  of  the  matrix  element  (wc  )"' . 

H.   LOWPASS  TO  HIGHPASS  ANALOG  TRANSFORMATION 

The  polynom.ial   substitution  for   the  lowpass  to  highpass 
transformation  is 


H(s)  =  Hp  (s) 


S  =  Wr  /S 


(2.31) 


Again      using      Equations      1.6      and      1.11,       the    first-order 

transf orm.ation   miatrix    is 


Ad)     = 


Wf 


1 
0 


(2.32) 


The  m:ethod   of  convolution  of  coefficients  can  be  used  to  get 


the  hiaher-order  transformation   matrices. 


However,  notice 


that  Equation  2.32  is  the  same  matrix  as  the  lowpass  to 
lowpass  transformation  with  the  reciprocal  taken  of  each  non 
zero  elemient,  and  then  each  row  reversed.  This  can  be  shown 
to  be  the  case  in  general,  so  that  the  lowpass  to  highpass 
transf  orm.ation  matrix  is  given  by 
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0 

0 


A(n)  = 


Wr  " 


0 

0 


1 1 


n  I 


(2.33) 


0 
0 


01 
I 

oi 
_1 


Again.  each  new  filter  coefficient  depends  only  on  one 
multiplication  although  the  calculation  of  the  matrix  element 
itself  will  require  calculations  involving  raising  the  design 
constant,  Wr  .  to  the  (n-1)^^  power  where  i  is  the  power  of  s 
in  the  final  filter  transfer  function. 

I.   LOWPASS  TO  BANDPASS  ANALOG  TRANSFORMATION 

The   polynomial   substitution   for  transforming  a  lowpass 
analog  prototype  to  a  bandpass  analog  filter  is 


H{s)  =  Hp  (s) 


s  = 


S2   +  W02 


Bs 


(2.34) 


where  wo  is  defined  as  the  geometric  mean  of  the  passband; 
Wo  =  (WuWi  ) '^  ,  and  B  is  the  width  of  the  passband,  B=Wu  -  wi  . 
(wu  and  wi  are  the  upper  and  lower  critical  frequencies 
respectively).   Again  using  Equations  1.6  and  1.11  and 
picking  out  the  appropriate  matrix  elements,  the  first-order 
transformation  matrix  is  given  by 
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B 


(2.35 


A(l)  = 


Wo 


ll 


The  convolution   of  coefficients   technique  is   again  used  to 

generate   the   higher   order  transformation  matrices  with  the 

results  for  the  first-order  to  third-order  matrices  listed  in 

Table  4.    As   with  the   bandpass  digital  transformation,  the 

analog  bandpass  transf orm.ation  will   double  the   order  of  the 

original  transfer   function  since  the  substitution  polynomial 

is  of  second  order. 

TABLE  4 
FIRST-ORDER  TO  THIRD-ORDER  BANDPASS 
ANALOG  TRANSFORMATION  MATRICES 


^ 


A'l)  =  I 


B 

0 


A(2)  =  10 


Wo 


B2 


0 


Bwo  2   0      B 
0      2wo  2   0 


01 

i 

0  i 

I 

I 
it 

_l 


A(3)  = 


0 

0 

i 
io 

I 

I 


B2  Wo  ^ 


Bwo  "   0 


B3 


B2 


2Bwo  2  0 


3v'r  "     0 


B 


3w-  2   0 


0 

0  i 
0  I 

i 

li 


Note  thc»t  each  r^atrix  is  slightly  over  half  empty  and 
that  the  nonzero  elements  are  sirple  function?  of  the  desian 
constants  Wo  and  B. 
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J.   LOWPASS  TO  BANDSTOP  ANALOG  TRANSFORMATION 

For  the  analog   bandstop   transformation,   the  polynomial 
substitution  is  the  reciprocal  of  the  bandpass  case: 


H(s)  =  Ho  (s) 


Bs 


<2.36) 


s  = 


S2   +  Wo 


Usina  the  same  method  as  before,  the  first-order 
transformation  matrix  is  found  to  be: 


Wo  2 

0 

1 

A(l)     = 

=   i 

0 

B 

oi 

l_ 

_J 

(2.37) 


Note  that  this  is  the  same  matrix  as  the  bandpass  case  but 
with  the  rows  in  reverse  order.   Using  the  convolution  of 
coefficients  technique  the  higher-order  matrices  are 
developed  and  the  results  for  the  first-order  to  third-order 
transf ormatiorj  matrices  are  shown  in  Table  5.  Notice  that  the 
reversal  of  the  order  of  the  rows  compared  with  the  bandpass 
case  is  true  m  general.   Therefore,  once  again  the  matrices 
are  over  half  empty  and  nonzero  elements  are  simple  functions 
of  B  and  Wo  . 
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TABLE  5 
FIRST-ORDER  TO  THIRD-ORDER  BANDSTOP 
ANALOG  TRANSFORMATION  MATRICES 


Wr 


A(l)     = 


B 


A{2)     = 


Wo 


2Wo2 


A(3)     = 


10 

1 

Bwo  " 

0 

B 

1 
0 

0 

B2 

0 

i 

Wo^ 

0 

2wo'' 

0 

0 

Bwo  " 

0 

2Bw 

0    • 

0 

B'  Wo  2 

0 

0 

0 

0 

B 

3wo2i; 


B2 


0 

1 

B 

0 

0 

0  i 

0 

1 

0 

1 

K.   GENERAL  OBSERVATIONS 

Now  that  the  various  transf orr.ation  matrices  for  digital 
filter  design  have  been  developed,  it  is  possible  to  compare 
the  two  methods  of  diaital  filter  design  mentioned  earlier 
and  depicted  in  Figure  1.  For  this  comparison,  consider  an 
analcc  lowpass  prototype  filter  transfer  function  with 
numerator  of  order  m,  and  denominator  of  order  n  (m  less  than 
or  equal  to  n) . 

First  consider  the  digital  to  digital  direct  design 
m.ethod.  Transf orm.ing  the  analog  prototype  to  a  digital 
prototvr-e  usina   the  bilinear  transf orm.ation  will  require  m+1 
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multiplications  and  additions  per  digital  prototype 
coefficient  for  the  numerator  polynomial,  and  n+1 
m^ultiplications  and  additions  per  digital  prototype 
coefficient  for  the  denominator  polynomial.  Then,  each  of 
the  digital  to  digital  transformations  will  require  2 (n+1) 
m.ultiplications  and  additions  since  a  result  of  the  bilinear 
transformation  is  to  make  n  =  m.  Therefore,  going  from  an 
analog  prototype  to  a  final  digital  filter  will  require 
(3n+m+4)  multiplications  and  additions.  If  lowpass  digital 
prototypes  are  developed  prior  to  the  design  process,  this 
can  be  reduced  to  2 (n+1)  since  then  only  the  digital  to 
digital  transformations  will  be  used  in  designing  the  final 
filter. 

Now  consider  the  analoa  design  method  for  the  same  case 
of  an  rr.'  ^  -order  numerator  and  n'  ^  -order  denominator  lowpass 
prototype  transfer  function.  For  the  lowpass  and  highpass 
case .  each  new  coefficient  will  require  only  1  multiplication 
and  addition.  For  the  bandpass  and  bandstop  transformations 
a  siirple  expression  for  the  number  of  calculations  per 
coefficient  is  difficult  to  find,  however  referring  to  Tables 
4  and  5,  the  number  of  multiplications  and  additions  per 
coefficient  is  always  less  than  half  the  order  of  the 
polynomial  to  be  transformed,  on  the  average.  This  gives 
approximately  (n-^^^.)/2  for  the  bandpass  to  bandstop  analog 
transf  orrr.ations .  Finally  the  analog  filter  must  be  converted 
to  a  dicital  filter   using  the   bilinear  transformation  which 
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will  require  (n+]T,+  2)  calculations  per  coefficient  for  the 
highpass  and  lowpass  case  and  2(n+F.+  l)  for  the  bandpass  and 
bandstop  case.  The  total  nunber  of  calculations  in  using  the 
analog  design  method  then  is  {n+K-t-4)  for  the  high/low  pass 
and  ( 2  .  5  (n+ir.) +2 )  for  the  bandpass/stop.  These  results  are 
sumriarized  in  Table  6. 


TABLE  6 
NUMBER  OF  MULTIPLICATIONS  AND  ADDITIONS  REQUIRED 
PER  COEFFICIENT  TO  FIND  THE  DIGITAL  FILTER 
COEFFICIENTS  FROM  A  LOWPASS  ANALOG  PROTOTYPE 
(WITH  AN  M^H  ORDER  NUMERATOR  AND  N^ «  ORDER  DENOMINATOR) 


Analog  Design 
high/low  pass    bandpass/stop 


analoc  transf orrr.at ions       '       2 

i 

.5 (n+m) 

i 

bilinear  transformation      1     n+m+2      1    2(n+m)+2      1 

1                      i                         ! 

i 

Total                   i     n+iii+4         2.5(n+ir.)+2 

1 

bilinear  transformation 


Digital  Design 
high/low  pass    bandpass/stop 


n+m+2 


n+iii+2m 


digital  transformations 


2  (n  +  l) 


2 (n+l) 


Total 


3(n  +  l)+m+l       I       3(n+l)+in+l 
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Referring  to  Table  6,  the  advantaae  of  one  method  over 
the  other  depends  upon  what  type  of  filter  is  being  designed 
and  the  order  of  the  prototype  transfer  functions  numerator 
and  denominator.  Consider  designing  a  fifth-order 
low/highpass  Chebyshev  or  Butterworth  filter-  In  this  case, 
the  order  of  the  numerator  polynomial  is  zero  and  the  analog 
method  will  have  a  slight  advantage  of  9  multiplications  and 
additions  versus  12  for  the  digital  design  with  precalculated 
prototypes  (and  19  for  digital  design  without  precalculated 
prototypes).  However,  if  a  bandpass /stop  elliptic  filter  was 
being  designed  with  n=5,  then  m=4  and  the  direct  digital 
design  method  (with  precalculated  prototypes)  has  an 
advantage  of  12  multiplications  and  additions  compared  to 
approximately  24  for  the  analog  method. 

A  final  consideration  in  favor  of  analog  design  is  the 
cor.plexity  of  the  elements  of  the  transformation  matrices. 
Since  the  elements  of  the  digital  transf orm^ation  matrices 
reguire  ir;uch  more  calculation  than  the  elements  of  the  analog 
transf  ormiation  matrices  or  the  bilinear  transformation 
(com.pare  Tables  1-5),  the  problem,  of  roundoff  error  is  much 
m.ore  serious  in  the  case  of  the  digital  to  digital  design 
method . 

L.   IMPLEMENTATION 

The  final  step  is  to  implement  the  algorithms  developed 
previously  in  a  computer  program;  to  assist  in  the  teaching  of 
digital  filter  design.   Since  the   program  is   intended  to  be 
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used  for  instructional  purposes,  the  nuiTierical  considerations 
discussed  earlier  are  not  considered  as  irr.portant  as 
conceptual  siir.plicity .  The  direct  design  method  was  chosen 
for  implementation  a?  the  design  constants  can  be  found 
easily  without  reference  to  "  prewarping  .'*  The  completed 
program.  called  'DFCADD'  (Digital  Filter  Computer  Aided 
Direct  Design) ,  is  contained  in  Appendix  A  with  its  own 
documentation.  The  program  follows  the  procedures  set  forth 
in  Chapter  10  of  Reference  2.  Prototype  coefficients  are 
stored  within  the  program:  for  first-order  to  fifth-order 
Butterworth  and  Chebyshev  (l/2,l,&2dB  ripples)  filters. 
While  second-order  to  sixth-order  elliptic  lowpass  prototype 
coefficients  are  calculated  within  the  program  for  any 
combination  of  1/2,1,2,  and  3  dB  passband  ripple  and  -20,- 
30 . -4C , -50 , -6C ,  and  -70  dB  stopband  attenuation,  with  the 
help  of  a  proararr,  developed  by  Professor  D.E.  Kirk  of  the 
Naval  Postgraduate  School  and  References  4  and  5.  As  a  final 
option,  the  user  m.ay  enter  his/her  own  prototype  coefficients 
for  up  to  a  tenth-order  analog  prototype.  After  calculating 
the  digital  filter  coefficients  an  option  exists  for  the 
rroara^  to  create  a  data  file,  titled  "filter",  with  the 
digital  filter  m.acnitude  and  phase  as  a  function  of  the 
diaital  frecuency  fror.  zero  to  pi.  This  data  file  can  then 
be  used  with  any  plottina  routine  to  generate  a  plot  of  the 
freauenv^v  resoonse  of  the  filter. 
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APPENDIX 
DFCADD  FORTRAN  PROGRAM 


FILE.  DFCADD    FORTRAN   A 


C 
C 

c 
c 
c 
c 
c 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 


9 
10 


X  THIS  PROGRAM  ASSISTS  IN  THE  CALCULATION  OF  DIGITAL  FILTER  X 

X  COEFFICIENTS  ACCORDING  TO  THE  PROCEDURE  SET  FORTH  IN  CHAPTER  x 

X  10  (SUMMARIZED  ON  PG  695)  OF  -  FIRST  PRINCIPLES  OF  DISCRETE  X 

X  SYSTEMS  AND  DIGITAL  SIGNAL  PROCESSING  -   BY  R.D.  STRUM   AND  » 

X  D.E.  KIRK  (REFERENCEC2)  )  X 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXKKXXXXXXXXXXX 

INTEGER  ORDER,TYPE, PASS, 0RDER2, END, AFORDF 

REAL  KPR,LAM,OMEGA(20),V(20),AO(20),BO(20),B1(20) 

COMPLEX  Z,NEWZ,ZN,ZD,HOFZN,HOFZD,HOFZ 

DOUBLE  PRECISION  TEMP,SUMR 

DOUBLE  PRECISION  R, SO, AOl , A02, BOl . B02, Bll, B12, HO 

DOUBLE  PRECISION  A03,  AO'i,  B03,  B0<4,  B13,  Bl<« 

DOUBLE  PRECISION  SFREQ, AFC, DFC, ALPHA, PI , PID^, DFCD2 

DOUBLE  PRECISION  ALC, AUC, DLC, DUC, SK, B, C, DLCD2, DUCD2 

DOUBLE  PRECISION  APROT, BILINR, DPROT, LPHPTR, BPBSTR, DFLTR 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXKXXXXXXXKX 

FORMAT  OF  ARRAYS; 

-APROT  CONTAINS  THE  ANALOG  PROTOTYPE 

APROT(NORD,PWR) 
-BILINR  IS  THE  ANALOG-DIGITAL  BILINEAR  TRANSFORMATION  MATRIX 

BILINR(ORDER,ROW,COL) 
-DPROT  IS  THE  DIGITAL  PROTOTYPE  FROM  APROT  AND  BILINR 

DPROT(N0RD,PWR) 
-LPHPTR  IS  THE  LP8HP  DIGITAL-DIGITAL  TRANFORMATION  MATRIX 

LPHPTR(ORDER,ROW,COL) 
-DFLTR  IS  THE  FINAL  DIGITAL  FILTER 

DFLTRCNORD,PWR) 
-BPBSTR  IS  THE  DP8BS  DIGITAL-DIGITAL  TRANSFORMATION  MATRIX 

BPBSTR(ORDER,ROW,COL) 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXX 

DIMENSION  APR0T(2,0:10) 
DIMENSION  BILINR(10,0:10,OilO) 
DIMENSION  DPR0T(2,0:10) 
DIMENSION  LPHPTRCIO, 0:10,0:10) 
DIMENSION  BPBSTRCIO, 0:10, 0:20) 
DIMENSION  DFLTR(2,0:20) 

XXXXXXXXXXXXXXXXXXXXXXKXXXXlOfXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

OPEN ( UNI T=8,FILE=' FILTER    DATA ', STATUS= 'UNKNOWN' ) 
xxxxxxxxxxxxxxxxxxxxxx*y»»>'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
PI  =  3.1'il5926535R98D+00 
SPI=3. 1^159 
PID'i  =  PI/<i.OD+00 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
CLEAR    APROT    PRIOR    TO    LOADING 
DO    1    K=l,2 

DO    2    L=0,10 
APROT(K,L)=0.0D+00 

CONTINUE 
CONTINUE 
XXXXXX^iXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

BEGIN  INTERACTION 

CALL  EXCMS  ( *CLRSCRN') 

HRITEt6,x)  'THIS  IS  A  PROGRAM  TO  CALCULATE  DIGITAL  FILTERS  USING' 

WRITE(6,x)  '  THE  DIRECT  DESIGN  METHOD.' 

HRITEC6,x) 

HRITE(6, X)   'XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX' 

WRITE(6,x)  'XHARNING-THIS  PROGRAM  IS  NOT  USER  FRIENDLY!!!  x' 
HRITE(6,x)  'x  YOU  MUST  ENTER  VALUES  OF  THE  APPROPRIATE  TYPEx' 
HRITEC6,x)  'X  (REAL  OR  INTEGER),  AND  PROPER  RANGE  X" 

l)RITE(6,X)   'XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX' 

HRITE(6,x) 
HRITE(6,x) 
READ(5,x) 


'   XX  ENTER  1  TO  CONTINUE  XX  • 

LOOK 

CALL  EXCMS  CCLRSCRN'  ) 

I1RITE(6,X)  'WHAT  TYPE  OF  FILTER  DO  YOU  WANTt' 
l)RITE(6,x)  '  (l-H  ALLOW  ONLY  UP  TO  5TH  ORDER  PROTOTYPES)' 
HRITE(6,x)  '       1-BUTTERllORTH  ' 

WRITE(6,x)  '       2-CHEBYSHEV  WITH  1/2  DB  RIPPLE  ' 
WRITE(6,x)  '       3-CHEBYSHEV  WITH  1  DB  RIPPLE  ' 
HRITE(6,x)  •      <i-CHEBYSHEV  WITH  2  DB  RIPPLE' 


NEWOOOIO 
NEW00020 
NEH00030 
NEW000<i0 
NEW00050 
NEH00060 
NEW00070 
NEM00080 
NEW00090 
NEHOOIOO 
NEHOQllO 
NEW00120 
NEI-I00130 
NEM00140 
NEW00150 
NEH00160 
NEW00170 
NEI'100180 
NEW00190 
NEH00200 
NEW00210 
NEW00220 
NEW00230 
NEW002^0 
NEH00250 
NEW0  026  0 
HEW00270 
NEH00280 
NEH00290 
NEW00300 
NEN00310 
NEW00320 
NEII00330 
NEH003<iO 
NEH00350 
NEW00360 
NEH00370 
NEII00380 
NEH00590 
NEHOOfiOO 
NEHOO'ilO 
NEW00«i20 
NEHOOaSO 
NEWOO^^O 
NEHOOISO 
NEH00<i60 
NEI100^70 
NEH00<<80 
NEN00^90 
NEN00500 
NEN00510 
NEW0a520 
NEII00550 
NEI1005<40 
NEII00550 
NEH00560 
NEH00570 
NEH00580 
NEN00590 
NEII00600 
NEH00610 
NEH00620 
NEH006  30 
NEW006^0 
NEH006  50 
NEH00660 
NEW00670 
NEI-100680 
NEN00690 
NEH00700 
NEH00710 
NEH00720 
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FILEi  DFCADD    FORTRAN   A 


5-ELLIPTIC  (ALLOWS  UP  TO  A  6TH  ORDER  PROTOTYPE) 
-PASSBAND  RIPPLE  0 . 5 , 1 . 0 , 2 . 0, OR  3.0  DB ' 
-STOPBAND  ATTENTUATION  -20, -30,-<i0, -50, -60,  • 
OR  -70  DB    • 
6-OTHER  -  ALLOWS  YOU  TO  USE  UP  TO  A  lOTH  ORDER' 
PROTOTYPE,  HOWEVER  YOU  MUST  ENTER  THE  ANALOG 
PROTOTYPE  COEFFICIENTS  • 

KX  INTEGER  1-6  «»• 
PE 

XKXKiOOOOCXXKXXXXXXXKXXXXXXKKKXKIIKKXKKXKMKKIIKXKKKKKXX 

•CLRSCRN') 

WHAT  IS  THE  ORDER  OF  THE  PROTOTYPE  YOU  HISH  TO  USET* 

-  MUST  BE  AN  INTEGER  IN  THE  RANGE  1-5  UNLESS  YOU' 
CHOSE  TO  ENTER  YOUR  OWN  PROTOTYPE  (THEN  1-10)' 
OR  YOU  CHOSE  AN  ELLIPTIC  FILTER  (THEN  1-6)' 

-  OR  ENTER  "0"  IF  YOU  WANT  HELP  SELECTING  THE  ORDER' 
BY  CHANGING  SPECIFICATION  FREQUENCIES  TO  THE' 
CORRESPONDING  DIGITAL  PROTOTYPE  FREQUENCIES  ' 

NOTEi  YOU  WILL  NEED  PROTOTYPE  FREQUENCY    ' 

RESPONSE  PLOTS  LIKE  THOSE  GIVEN  IN  REF(2)  FOR   ' 

BUTTERHORTH  FILTERS  (FIG  10.33,  PGS  6918692)   • 

AND  CHEBYSHEV  FILTERS  (FIG  10.3<4,  PGS  696-698)  • 

DER 

^). AND. (TYPE. LT.8))THEN 
DTHEN 
•ELLIPTIC  FILTERS  GIVEN  FOR  2ND-9TH  ORDER  ONLY' 


WRITE{6,») 
WRITE(6,«) 
WRITE(6,*) 
WRITE(6,x) 
WRITE(6,K) 
WRITE(6,K) 
WRITE(6,x) 
WRITE(6,*) 
WRITE(6,K) 
READ(5,K)  TY 

C       XXXXXXKXXXXX 

l<i    CALL  EXCMS  ( 
15    WRITE(6,x) 
WRITE(6,x) 
WRITE(6,x) 
WRITE(6,x) 
WRITE(6,X) 
WRITE(6,X) 
WRITE(6,K) 
WRITE(6,x) 
WRITE(6,K) 
WRITE(6,x) 
HRITE(6,X) 
HRITE(6,x) 
HRITE(6,*) 
READ(5,x)  OR 
IF((TYPE.GT 
IF(ORDER.EQ 
WRITE(6,K) 
WRITE(6,x) 
GO  TO  10 
ENDIF 
ENDIF 
C 

CALL  EXCMS  ('CLRSCRN') 

HRITE(6,x)  'DO  YOU  WANT  TO  ENTER  SPECIFICATION  FREQUENCIES' 

HRITE(6,x)  'AS  ANALOG  OR  DIGITAL  FREQUENCIES?' 

HRITE(6,X)  '     1-ANALOG   ' 

HRITE(6,x)  '     2-DIGITAL  ' 

HRITE(6,x)  '   XX  INTEGER  1  OR  2  xx  • 

READ(5,x)  AFORDF 

CALL  EXCMS  ( 'CLRSCRN'  ) 

IF(AFORDF.EQ. DTHEN 

HRITE(6,X)  'WHAT  IS  THE  SAMPLING  FREQUENCY  IN  KHZ?' 
NRITE(6,x)  '    xx  DECIMAL  xx' 
READ(5,x)  SSFREQ 
SFREQ=DBLE(SSFREQ) 
ENDIF 

CALL  EXCMS  ( 'CLRSCRN') 
HRITE(6,X) 
11RITE(6,X) 
HRITE(6,x) 
HRITE(6,X) 
HRITE(6,X) 
IIRITE(6,X) 
HRITE(6,X) 
HRITE(6,X)  ' 
READ(5,x)  PASS 
C     FORMULAS  FOR  HIGHPASS  AND  LOWPASS 
IF((PASS.E(3.1)  .0R.(PASS.EQ.2))THEN 
CALL  EXCMS  ('CLRSCRN') 
IFCAFORDF.EQ. DTHEN 

IIRITE(6,X)  'WHAT  IS  THE  CRITICAL  FREQUENCY  IN  KHZ?  ' 
IIRITE(6,x)  '     -  NEG  3DB  PT  FOR  BUTTERWORTH   • 
HRITE(6,x)  '     -  RIPPLE  EDGE  FOR  ELLIPTIC  AND  CHEBYSHEV 
WRITE(6,x)  '        XX  DECIMAL  xx  • 
READ(5,x)  SAFC 
AFC=DBLE(SAFC) 

DFC=(AFC/SFREQ)XPIX(2.0D+00) 
ENDIF 
1F(AF0RDF.EQ.2)THEN 

WRITE(6,x)  'WHAT  IS  THE  CRITICAL  DIGITAL  FREQUENCY?' 


WHAT  TYPE  OF  PASSBAND  DUE  YOU  WANT?' 
1-LOWPASS' 
2-HIGHPASS' 
3-BANCPASS' 
<i-BANDSTOP' 

XX  INTEGER  l-<i  xx  • 


•NEW00730 
NEW007<40 
NEW00750 
NEW00760 
NEW00770 

'NEW00780 
NEW00790 
NEW00800 
NEW00810 
NEW00820 
NEW00830 
NEW008<40 
NEW00850 
NEW00860 
NEW00870 
NEW00880 
NEW00890 
NEW00900 
NEW00910 
NEW00920 
NEW00930 
NEW009^0 
NEW00950 
NEW00960 
NEW00970 
NEW00980 
NEW00990 
NEWOIOOO 
NEWOIOIO 
NEW01020 
NEW01030 
NEW010<<0 
NEW01050 
NEH0106  0 
NEW0107  0 
NEM01080 
NEW01090 
NEHOllOO 
NEWOlllO 
NEH01120 
NEW01130 
NEII011<iO 
NEllonSO 
NEW01160 
NEW0117  0 
NEII01180 
NEII01190 
NEH01200 
NEII01210 
NEW01220 
NEH01230 
NEW012<iO 
NEII01250 
NEII01260 
NEW0127  0 
NEW01280 
NEW01290 
NEN01300 
NEW01310 
NEH01320 
NEH01330 
NEW013<iQ 
NEW01350 
NEW01360 
NEMO  137  0 
NEW01380 
NEW01390 
NEW01'<00 
NEW01<il0 
NEW01<<20 
NEW01<i30 
NEW01<>^0 
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FILE.  DFCADD 


FORTRAN   A 


C 
C 


c 

c 


c 
c 


HRITE(6,K)  •      -  NEG  3DB  FOR  BUTTERWORTH  ' 
WRITE(6,K)  •      -  RIPPLE  EDGE  FOR  ELLIPTIC  AND  CHEBYSHEV 
HRITE(6,»)  '       «K  DECIMAL  KK  • 
READ(5,X)  SDFC 
DFC=DBLE(SDFC) 
ENDIF 

DFCD2=DFC/2.0D+00 
CALCULATE  THE  APPROPRIATE  VALUES  OF  THE  DESIGN  CONSTANT  ALPHA 
FROM  TABLE  10.7  REFERENCE  (2)  AND  REFERENCE  (3) 
IFCPASS.EQ.DTHEN 

ALPHA  =  DSIN((PID<i)-(DFCD2))/DSIN((PID'i)  +  (DFCD2)) 
ENDIF 
IF(PASS.EQ.2)THEN 

ALPHA  =  DC0S((PID«i)-(DFCD2))/DC0S((PID<i)  +  (DFCD2)) 
ENDIF 

OPTION  TO  SHOW  CALCULATED  VALUE  OF  ALPHA 
CALL  EXCMS  CCLRSCRN') 

HRITE<6,X)  -DO  YOU  WISH  TO  SEE  THE  CALCULATED  VALUE  OF  ALPHA?' 
WRITE(6,X)  •  XX  1-YES/2-N0  xx' 
READ(5,x)  LOOK 
IFCLOOK.EQ.DTHEN 

WRITE(6,*) 

WRITE(6,903)  ALPHA 

1IRITE(6,») 

WRITE(6,*)  'XX  ENTER  1  TO  CONTINUE  xx • 

READ(5,«)  LOOK 
ENDIF 
ENDIF 

XXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
BANDPASS  AMD  BANDSTOP  FORMULAS 
IF((PASS.EQ.3)  .OR.(PASS.EQ.<f))THEN 
CALL  EXCMS  (  'CLRSCRN') 
IFCAFORDF.EO.DTHEN 

WRITE(6,x)  'WHAT  IS  THE  LOWER  CRITICAL  FREQUENCY  IN  KHZT' 
'    -  NEG  3DB  FOR  BUTTERWORTH' 

'    -  RIPPLE  EDGE  FOR  ELLIPTIC  AND  CHEBYSHEV   ' 
'       XX  DECIMAL  XX' 
SALC 


WRITE(6,X) 

WRITE(6,x) 

WRITE(6,X) 

READ(5,K) 

ALC=DBLE(SALC) 

CALL  EXCMS  (  'CLRSCRN') 


'WHAT  IS  THE  UPPER  CRITICAL  FREQUENCY  IN  KHZ?' 

'    -  NEG  3DB  FOR  BUTTERWORTH' 

'    -  RIPPLE  EDGE  FOR  ELLIPTIC  AND  CHEBYSHEV  ' 
'       XX  DECIMAL  XX' 


WRITE(6,*) 
WRITE(6,X) 
WRITE(6,x) 
WRITE(6,x) 
READ(5,x)  SAUC 
AUC=DBLE(SAUC) 

DLC=(ALC/SFRE0)»(2.0D+0  0)»PI 
DUC=(AUC/SFREQ)X(2.0D+0  0JXP1 
ENDIF 
IF(AF0RDF.E0.2)THEN 

•WHAT  IS  THE  LOWER  CRITICAL  DIGITAL  FREQUENCY?' 
'    -  NEG  3DB  FOR  BUTTERWORTH  ' 
'    -  RIPPLE  EDGE  FOR  CHEBYSHEV  AND  ELLIPTIC 
XX  DECIMAL  XX  ' 


WRITE(6,x) 

WRITE(6,*) 

WRITE(6,X) 

WRITE(6,X) 

READ(5,*)  SDLC 

DLC=DBLE(SDLC) 

CALL  EXCMS  CCLRSCRN') 


WHAT  IS  THE  UPPER  CRITICAL  DIGITAL  FREQUENCY?' 

-  NEG  3DB  FOR  BUTTERWORTH  ' 

-  RIPPLE  EDGE  FOR  CHEBYSHEV  AND  ELLIPTIC 
XX  DECIMAL  XX  • 


WRITE(6,x) 
WRITE(6,x) 
WRITE(6,x) 
WRITEr6,x) 
READ(5,x)  SDUC 
DUC=DBLE(SDUC) 
ENDIF 

DLCD2=DLC/2.0D+00 
DUCD2=DUC/2.0D+00 
CALCULATE  APPROPRIATE  VALUES  OF  DESIGN  CONSTANTS  ALPHA  AND  K 
FROM  TABLE  10.7  REFERENCE  (2)  AND  REFERENCE  (3) 
ALPHA=DC0S((DUCD2)+(DLCD2))/DC0S((DUCD2)-(DLCD2)) 
IF(PASS.EO.«)THEII 
RK  =  DTAN(PID'i)XDTAN((DUCD2)-(DLCD2)) 


NEWOl'iSO 
NEW01<<60 
NEW01<»70 
NEW01<i80 
NEW01<490 
NEW01500 
NEW01510 
NEW01520 
NEW01530 
NEWOlS'iO 
NEW01550 
NEW01560 
NEW01570 
NEW01580 
NEW01590 
NEW01600 
NEW01610 
NEW01620 
NEW01630 
NEW016<*0 
NEW01650 
NEW01660 
NEW01670 
NEW01680 
NEW01690 
NEW01700 
NEW01710 
NEW01720 
NEW01730 
NEWOW'iO 
NEW01750 
NEW01760 
NEW01770 
NEW01780 
NEW01790 
NEH018  00 
NEW01810 
NEW01820 
NEW01830 
NEWOIB^O 
NEW01850 
NEW0186  0 
NEW0187  0 
NEW01880 
NEW01890 
NEW01900 
NEW01910 
NEW01920 
NEH01930 
NEWOn^iO 
NEW01950 
NEW01960 
NEW01970 
NEW0198  0 
NEW01990 
NEW0200  0 
NEH02010 
MEH02020 
HEW02030 
NEW020'iO 
NEW02050 
NEW02060 
NEW02070 
NEW02080 
NEW02090 
NEW02100 
NEW02110 
NEW02120 
NEW02130 
NEW021<iO 
NEW02150 
NEW02160 
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B=((2.0D+00)KALPHA)/(RK+1.0D+00)  NEW0217  0 

C=((1.0D+00)-RK)/((1.0D+00)+RK)  NEW02180 

ENDIF  NEH02190 

1F(PASS.EQ.3)THEN  NEW02200 

RK  =  DTAN(PID<i)/DTAN((DUCD2)-(DLCD2))  NEW02210 

B=((2.0D+00)KALPHA»RK)/(RK+(1.0D+00))  NEW02220 

C=(RK-(1.0D+00))/(RK+(1.0D+00))  NEW02230 

ENDIF  NEW022<iO 

C     OPTION  TO  SHOW  CALCULATED  VALUE  OF  ALPHA  AND  K  NEW02250 

CALL  EXCMS  CCLRSCRN')  NEW02260 

WRITE(6,»)  'DO  YOU  WISH  TO  SEE  THE  CALCULATED  VALUE  OF  ALPHA  •    NEW02270 

WRITE(6,»()  '   AND  KT   •  NEW02280 

WRITE(6,x)  •      x«  1-YES/2-N0  ««•  NEW02290 

READ(5,»)  LOOK  NEW02300 

IFCLOOK.EQ.DTHEN  NEW02310 

WRITE(6,K)  NEW02320 

WRITE(6,903)  ALPHA  NEW02330 

WRITE(6,90<i)  RK  NEH023<iO 

WRITE(6,X)  NEW02350 

WRITE(6,»)  'XX  ENTER  1  TO  CONTINUE  XX'  NEII02360 

READ(5,X)  LOOK  NEW02370 

ENDIF  NEW02380 

ENDIF  NEW02390 

C       XXKXXXXXXXXXXXXXXXXXXKXXKXXXKKXXXXXKXXXXXXXXXXXXXXXKXXXKXXXXXXXX   NEH02^00 

C     OPTION  TO  ASSIST  IN  DETERMINING  THE  PROTOTYPE  ORDER  BY  CHANGING   NEW02'ilO 

C     SPECIFICATION  FREQUENCIES  TO  APPROPRIATE  PROTOTYPE  FREQUENCIES    NEW02^20 

C     FROM  EQN  10.27^  REFERENCECl)  NEW02<i30 

IF(ORDER.EQ.0)THEN  NEW02'i'iO 

16     CALL  EXCMS  CCLRSCRN')  NEW02'i50 

IFCAFORDF.EQ.DTHEN  NEH02'460 

WRITE(6,x)  'WHAT  IS  THE  ANALOG  SPECIFICATION  FREQUENCY  (KHZ)'   NEH02'i70 

HRITE(6,x)  'THAT  YOU  WANT  TO  CONVERT  TO  A  DIGITAL*  NEW02^80 

WRITE(6,x)  'LOWPASS  PROTOTYPE  FREQUENCY?'  NEW02'i90 

WRITE(6,*)  '   XX  DECIMAL  xx •  NEN02500 

READ(5,x)  F  NE1102510 

DF=(2.0)XSPIX(F/SSFREQ)  NEH02520 

ENDIF  HEH02530 

IF(AF0RDF.EQ.2)THEN  NEl)025'40 

NRITE(6,X)  'WHAT  IS  THE  DIGITAL  SPECIFICATION  FREQUENCY'  NEW02550 

HRITE(6,x)  'THAT  YOU  WANT  TO  CONVERT  TO  A  DIGITAL'  NEW02560 

WRITEt6,x)  'LOWPASS  PROTOTYPE  FREQUENCY?'  NEW02570 

WRITE(6,x)  •   XX  DECIMAL  xx'  NEW02580 

READC5,x)  DF  NEW02590 

ENDIF  NEW02600 

X  =  COS(DF)  NEII02610 

Y=SIN(DF)  NEW02620 

Z=CMPLX(X,Y)  NEW02630 

IF(PASS.EQ.1)THEN  NEW026^0 

AS=REAL(ALPHA)  NEW02650 

NEIIZ=(Z-AS)/((1 .O)-(ASXZ))  NEW02660 

ENDIF  NEW02670 

IF(PASS.EQ.2)THEN  NEW02680 

AS=REAL(ALPHA)  NEW02690 

NEWZ=(-1.0)X((Z-AS)/(C1.0)-(ASXZ)))  NEW027  00 

ENDIF  NEW02710 

IF(PASS.EQ.3)THEN  NEW02720 

BS=REAL(B)  NEW02730 

CS  =  REAL(C)  NEW027<iO 

NEIIZ=(-1.0)X(((ZXX2)-(BSXZ)+CS)/(C1.0)-(BSXZ)+(CSX(ZXX2))))  NEW027  50 

ENDIF  NEW02760 

IF(PASS.EQ.<i)THEN  NEW02770 

BS=REAL(B)  NEW02780 

CS=REAL(C)  NEW02790 

NEWZ=((ZXK2)-(BSXZ)+CS)/((1.0)-(BSXZ)+(CSX(ZXX2)))  NEW028  00 

ENDIF  NEW02810 

PTHETA=ABS(ATAN(AIMAGCNEWZ)/REAL(NEWZ)))  NEW02820 

IFCPASS.EQ.DTHEN  NEW02830 

IFCF.GT.SAFOTHEN  NEW028<iO 

PTHETA=SPI-PTHETA  NEW02850 

ENDIF  NEW02860 

ENDIF  NEW02870 

IF(PASS.EQ.2)THEN  NEW02880 
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IFCF.LT.SAFOTHEN  NEW02890 

PTHETA=SPI-PTHETA  NEW02900 

ENDIF  NEH02910 

EHDIF  NEW02920 

IF(PAS5.EQ.3)THEN  NEW02950 

IF((F.LT.SALC).OR.(F.GT.SAUC))THEN  NEW029<»0 

PTHETA=SPI-PTHETA  NEW02950 

ENDIF  NEW02960 

ENDIF  NEW02970 

IF(PASS.E(J.<i)THEN  NEW02980 

IF((F.GT.SALC).AND.(F.LT.SAUC))THEN  NEW02990 

PTHETA=SP1-PTHETA  NEW03000 

ENDIF  NEH03010 

ENDIF  NEW03020 

WRITE(6,905)  PTHETA  NEW03030 

WRITE(6,x)  NEW030<iO 

WRITE(6,x)  'WANT  TO  CONVERT  ANOTHER  FREQUENCYt '  NEH03050 

WRITE(6,x)  •   «x  1-YES/2-N0  xx'  NEW03060 

READ(5,x)  LOOK  NEW03070 

IFCLOOK.EQ.DGO  TO  16  NEW03080 

IF(LOOK.E(5.2)THEN  NEW03090 

HRITE(6,x)  •  NEW03100 

WRITE(6,x)  'WHAT  ORDER  PROTOTYPE  DO  YOU  WANT  TO  USE?'  NEW03nO 

WRITE(6,x)  '   XX  INTEGER  1-5  xx  •  NEW03120 

READ(5,x)  ORDER  NEW03I30 

ENDIF  NEW03I«(0 

ENDIF  NEW03150 

C       XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXKXXXXXX   NEW03I6  0 

C      LOAD  APROT  MATRIX  WITH  APPROPRIATE  COEFFICIENTS  NEW03170 

C      APROT(TYPE=I)  ARE  BUTTERWORTH  FILTERS  NEW03180 

IFCTYPE.EQ.DTHEN  NEW03190 

IFCORDER.EQ.DTHEN  NEII03200 

APROT(1,0)=1.0D+00  NEW03210 

APROT(2,0)=1.0D+00  NEW03220 

APROT(2,I)=1.0D+00  NEW03230 

ELSEIF(ORDER.E(5.2)THEN  NEH032'40 

APROT(1,0)=I.OD+00  NEW03250 

APROT(2,0)=I .OD+00  NEW03260 

APROT(2a)  =  D$ORT(2.0D+00)  NEW0327  0 

APROTC2,2)=1.0D+00  NEW03280 

ELSEIF(0RDER.E0.3)THEN  NEW03290 

APROT(I,0)=1 .OD+00  NEH03300 

APROT(2,0)=1.0D+00  NEN033I0 

APR0T(2,I)=2. OD+00  NEW03320 

APR0T(2, 2)^2. OD+00  NEW03330 

APR0T(2,3)=1. OD+00  HEH033^0 

ELSEIF(ORDER.EQ.'i)THEN  NEW03350 

APROT(I,0)=1. OD+00  NEW03360 

APRaT(2, 0)^1 .OD+00  HEW03370 

APROT(2,1)=2.613126D+00  HEW0338  0 

APR0T(2,2)  =  3.<H^21<iD+0  0  MEW03390 

APROT(2,3)=2.613126D+00  NEW03'i00 

APR0T(2,<i)  =  l  .OD+00  NEW03<ilO 

ELSE1F(0RDER.EQ.5)THEN  NEH03^20 

APR0T(1,0)  =  1  .OD+00  NEWOSiiSO 

APROT(2,0)  =  1. OD+00  NEHOS-i^O 

APROT (2, 1)  =  3.236  068  0+00  NEH0  3<i50 

APROT  (2,  2)  =5.236  0680+ 00  NEII03^6  0 

APROT (2,  3)  =5.236  068  0+ 00  NEH03'^7  0 

APR0T(  2,  li)  =  3.236  0680+00  NEW03'i80 

APR0T(2,5)  =  1. OD+00  NEII03^90 

EHDIF  NEH03500 

ENDIF  NEW035I0 

C      APR0T(TYPE=2)  ARE  CHEBYSHEV  FILTERS  WITH  .5DB  RIPPLE  NEII03520 

IF(TYPE.EQ.2)THEN  NEW03530 

IFCORDER.EQ.DTHEN  NEW035<iO 

APROT (1,0) =2.8627750+00  NEW03550 

APROT (2,0) =2. 86 277 50+00  NEW03560 

APR0T(2,1)=1 .00+00  NEW03570 

ELSEIF(0RDER.EQ.2)THEN  HEH03580 

APR0T(l,0)  =  I.<i31388D+00  NEH0  3590 

APROT(2,0)=I.516205D+00  NEW03600 
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.379051D+00 
.025<i55D+00 
.716866D+00 
.197386D+00 
•OD+OO 


.178923D+00 
,7525180+00 
.3095750+00 
.9373680+00 
.172'i91D+00 
.00+00 


APROT(2,l)  =  l.<i25625D+00 

APR0T(2,2)=1. 00+00 
ELSEIF(0RDER.EQ.3)THEN 

APROT(l,0)  =  0.71569<iD+00 

APROT(2,0)  =  0.71569<iD+00 

APROT( 2,1) =1.53^8950+00 

APROTC 2, 2)  =  1.2529130+00 

APR0T(2,3)=1. OD+OO 
ELSEIF(OROER.EQ.<i)THEN 

APR0T(l,0)=O.3578'i7O+00 

APROT(2,0)=0. 

APR0T(2,1)=1 

APR0T(2,2)=1 

APR0T(2,3)=1 

APR0T(2,<i)  =  l 
ELSEIF(0R0ER.E(}.5)THEN 

APROT(1,0)=0. 1789230+00 

APROT(2,0)=0 

APROT(2,1)=0, 

APR0T(2,2)=1. 

APR0T(2,3)=1, 

APR0T(2,'1)  =  1 

APR0T(2,5)=1, 
ENDIF 
ENDIF 

APR0T(TYPE=3)    ARE   CHEBYSHEV 
1F(TYPE.EQ.3)THEN 
IFCORDER.EO.DTHEN 

APROT(1,0)=1. 1965230+00 

APROT(2,0)=1. 1965230+00 

APR0T(2,1)=1 .00+00 
ELSEIF(0RDER.EQ.2)THEN 

APROT(1,0)=0.982613D+00 

APROT(2,0)=1 .1025100+00 

APROT(2,l)  =  1.0977  3'iO+0  0 

APR0T(2,2)=1. 00+00 
ELSEIF(0RDER.EQ.3)THEN 

APR0T(1,0)  =  0. '4913070+00 

APR0T(2,0)=0. 4913070+00 

APROK  2,  l)  =  1.238'i  090+00 

APR0T(2,2)  =  0.9883<;iO+00 

APR0T(2,3)=1 .00+00 
ELSEIF(ORDER.EQ.'i)THEN 

APROT(1,0)=0. 2456530+00 

APROT(2,0)=0. 

APR0T(2,1)=0 

APR0T(2,2)=1 

APR0T(2,3)=0 

APR0T(2,4)=1 
ELSEIFCORDER. 


.2756280+00 
.7426190+00 
.4539250+00 
.9528110+00 
.00+00 
.EQ.5)THEN 


FILTERS  WITH  lOB  RIPPLE 


APROTCl, 0=0.1228270  +  00 

APROT(2,0) =0.122827  0+00 

APROT(2,1)=0. 5805340+00 

APROT(2,2)=0. 9743960+00 

APROK 2, 3) =1.688816 0+00 

APROT(2,4)=0.936  8  20D+00 

APR0T(2,5)=1. 00+00 
ENDIF 
ENDIF 

APR0T(TYPE=4)    ARE   CHEBYSHEV    FILTERS   WITH    2DB    RIPPLE 
IF(TYPE.EQ.4)THEN 
IFCOROER.EO.DTHEN 

APROTCl, 0)=1. 3075600+00 

APROT(2,0)=1. 3075600+00 

APR0T(2,1)=1. OD+OO 
ELSEIF(0R0ER.EQ.2)THEN 

APROK 1,0) =0.5058 030+00 

APROT(2,0)=0.636768D+00 

APROT(2,1)=0. 8038160+00 

APR0TC2,2)=1. 00+00 
ELSE1F(0RDER.EQ.3)THEN 

APROT(1,0)=0. 3268900+00 


NEW03610 
NEI'103620 
NEW03630 
NEW03640 
NEW03650 
NEM03660 
NEW03670 
NEW03680 
NEN03690 
NEW03700 
NEW03710 
NEW03720 
NEW03730 
NEW03740 
NEW03750 
NEW03760 
NEW03770 
NEW03780 
NEH0379D 
NEW03800 
NEW03810 
NEW03820 
NEW03830 
NEt'103840 
NEW03850 
NEtl03860 
NEW03870 
NEH03880 
NEH03890 
NEH03900 
NEH03910 
NEH0392a 
NEH03930 
NEtl03940 
NEW03950 
NEH03960 
NEH03970 
NEII03980 
NEH03990 
NEII04000 
NEH04010 
NEH04020 
NEH04030 
NEI'104040 
NEH04050 
NEH04060 
NEH04070 
NEH04080 
NEH04090 
NEH04100 
NEH04110 
NEH04120 
NEII04130 
NE1I04140 
NEH04150 
NEI'104160 
NEW04170 
NEH04180 
NEtl04190 
NEH04200 
NEII04210 
NEH04220 
NEW04230 
NEH04240 
NEH04250 
NEW04260 
NEH04270 
NEI-104280 
NEH04290 
NEII04300 
NEH04310 
NEW04320 
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C 
C 

c 
c 


APROTC 2,0) =0.3268900+ 00 

APROT( 2,1) =1.0221900+00 

APROTC 2, 2) =0.7 378220+00 

APR0T(2,3)=1. 00+00 
ELSEIF(ORDER.EQ.<i)THEN 

APROT(1,0)  =  0.16  3<*^50+00 

APROTC 2,0) =0.20576 50+00 

APR0T(2,1) =0.5167980+00 

APROT(2,2)  =  1.256'i82D+00 

APROTC 2, 3) =0.7162150+00 

APR0TC2,^)=1. 00+00 
ELSEIFC0RDER.EQ.5)THEN 

APROTC 1,0) =0.08 17230+00 

APROTC2,0)=0. 0817230+00  ' 

APROTC2,l)  =  0.<i593^90+00 

APROTC 2, 2) =0.6 93477 0+00 

APROTC2,3)  =  l.'i995430+00 

APROTC 2, 4) =0.7 06 46 10+00 

APR0TC2,5)=1. 00+00 
ENDIF 
ENOIF 

APR0TCTYPE=5)  ARE  ELLIPTIC  FILTERS 

ELLIPTIC  COEFFICIENTS  ARE  CALCULATEO  USING  A  PROGRAM  DEVELOPED 
BY  PROFESSOR  0.  E.  KIRK,  NAVAL  POSTGRADUATE  SCHOOL,  MONTEREY, 
CALIFORNIA   93943 
A02=0. 00+00 
B02=0. 00+00 
612=0.00+00 
S0=0. 00+00 
IFCTYPE.E0.5)THEN 

CALL  EXCMS  C  'CLRSCRN') 

HRITEC6,X)  'WHAT  PASSBAND  RIPPLE  DO  YOU  WANT?  CIN  OB)' 

WRITEC6,X)  'xxx  MUST  BE  EITHER  0.5,  1.0,  2.0,  OR  3.0  »x»' 

REA0(5,«)  ASUBP 

HRITEC6,x) 


•WHAT  STOPBANO  ATTENUATION  00  YOU  WANT?  CIN  OB)' 
•XXX  MUST  BE  -20,-30,-40,-50,-60,  OR  -70  «xx' 
•XXX  00  NOT  INCLUDE  A  DECIMAL  xxx' 


WRITE(6,X) 
WRITE(6,x) 
WRITEC6,x) 
WRITEC6,x) 
REA0C5,x)  NATTN 
NATTN=-1XNATTN 
ASUBA=REALCNATTN) 
IFCASUBP.E0.0.5)THEN 
IFCNATTN.EQ.20)THEN 

IFC0RDER.EQ.2)SR=2. 76261 
IFC0RDER.EQ.3)SR=1 .42189 
IFC0RDER.EQ.4)SR=1 .15188 
IFC0RDER.EQ.5)SR=1 .04465 
I FC ORDER. EQ. 6 )5R= 1.01553 
ENOIF 
IFCNATTN.E0.30)THEN 

IFC0RDER.EQ.2)SR=4. 80880 
IFC0RDER.EQ.3)SR=1 .92322 
IFCaRDER.EQ.4)SR=l. 32446 
IFC0RDER.EQ.5)SR=1. 12912 
IFCORDER.EQ.6)SR=1.05394 
ENOIF 
IFCNATTN,EQ.40)THEN 

I FC ORDER. EQ. 2 )SR=8. 848925 
IFC0R0ER.EQ.3)SR=2. 71147 
IF(0RDER.E0.4)SR=1 .62842 
IFC0RDER.EQ.5)SR=1. 27264 
IFCOROER. Eg. 6)SR=1. 12697 
ENOIF 
IFCNATTN.EQ.50)THEN 

IFC0RDER.EQ.2)SR=15.06  069 
IFCORDER.EQ.3)SR=3.9  0430 
IFCOROER. EQ.4)5R=2. 06924 
IFC0R0ER.EQ.5)SR=1.4846  9 
IFCOROER.E(3.6)SR  =  1.24101 
ENDIF 
IFCNATTN.EQ.60)THEN 


NEW04330 
NEW04340 
NEW0435a 
NEW04360 
NEW04370 
NEW04380 
NEW04390 
NEW04400 
NEW04410 
NEW04420 
NEW04430 
NEW04440 
NEW04450 
NEW04460 
NEW04470 
NEW04480 
NEW04490 
NEW04500 
NEW04510 
NEW04520 
NEW04530 
NEW04540 
NEW04550 
NEW04560 
NEW04570 
NEW04580 
NEW04590 
NEW04600 
NEW04610 
NEW04620 
NEW04630 
NEW04640 
NEW04650 
NEW04660 
NEW04670 
NEH04680 
NEW04690 
NEW04700 
NEW04710 
NEH04720 
NEW04730 
NEW04740 
NEW04750 
HEW04760 
NEW04770 
NEH04780 
NEW04790 
NEII04800 
NEW04810 
HEW04820 
NEW04830 
NEH04840 
NEW04850 
NEW04860 
NEW0487  0 
NEW04880 
NEH04890 
NEW04900 
NEN04910 
NEW04920 
NEW04930 
NEW04940 
NEW04950 
NEW04960 
NEH04970 
NEN04980 
NEH04990 
NEW05000 
NEW05010 
NEH05020 
NEN05030 
NEW05040 
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ENDIF 

IF(NATTN.E(5.<iO)THEN 

IF(0RDER.EQ.2)SR=5.76107 
IF( ORDER. EQ. 3 )SR=2. 13923 
IF(  ORDER.  EQ.<«)SR  =  l.'(08<i2 
IF(0RDER.EQ.5)SR=1. 16811 
IF(ORDER.EQ.6)SR=1.07316 

ENDIF 

IF(NATTN.E(}.50)THEN 

IF(ORDER.E().2)SR  =  10.19181 
IF(ORDER.E0.3)SR=3.0^137 
IF(0RDER.EQ.<i)SR  =  1.7  5285 
IF(0RDER.EQ.5)SR=1.332^3 
IF(0RDER.EQ.6)SR=1.15868 

ENDIF 

IF(NATTN.EQ.60)THEN 

IF(ORDER.E(5.2)SR  =  18.09398 
IF(0RDER.EQ.3)SR  =  <».39729 
IF(ORDER.EQ.^)SR  =  2.2'('i^0 
IF(0RDER.EQ.5)SR=1. 56860 
IF(GRDER.EQ.6)SR=1.28693 

ENDIF 

IF(NATTN.EQ.70)THEN 

IF(0RDER.EQ.2)SR=32.15951 
IF(ORDER.EQ.3)SR  =  6.<i089^ 
IF(0RDER.EQ.^)SR=2.9236  3 
IF(0RDER.EQ.5)SR=1 .889  06 
IF(ORDER.EQ.6)SR  =  l.'i6330 

ENDIF 
ENDIF 
IF(ASUBP.E0.3.0)THEN 

IFCNATTN.EQ.20)THEN 

IF(0RD£R.EQ.2)SR=1.7  3915 
IF(0RDER.EQ.3)SR=1 .15516 
I F( ORDER. EQ.^)5R= 1.03853 
IF(0RDER.EQ.5)SR=1. 00996 
IFCORDER. Eg. 6)SR=1. 00260 

ENDIF 

IF(NATTN.EQ.30)THEN  ' 
IF(0RDER.EQ.2)SR=2. 90352 
IFtORDER.EQ.3)SR  =  l.<i581<i 
IF(0RDER.EQ.'i)SR  =  l.l<i5a2 
IF(ORDER.EQ.5)SR=l.n5020 
IF(ORDER.EQ.6)SR5l.01783 

ENDIF 

IF(NATTN.EQ.^O)THEN 

IF(ORDER.EQ.2)SR=5.0558^ 
IF(0RDER.EQ.3)SR=1. 98022 
IF(0RDER.EQ.^)SR  =  1.3'i663 
IF(0RDER.EC!.5)SR  =  1  .1393<i 
1F(0RDER.EQ.6)SR=1. 05892 

ENDIF 

IF(NATTN.EQ.50)THEN 

IF(0RDER.EQ.2)SR=8. 93009 
IF(0RDER.EQ.3)SR=2. 79862 
IF(0RDER.EQ.'i)SR  =  1.661'48 
IF(0RDER.E5.5)SR=1. 28850 
IF(0RDER.EQ.6)SR  =  1.1353<i 

ENDIF 

IF(NATTN.EQ.60)THEN 

IF(0RDER.EQ.2)SR=15. 8^610 
IF(0RDER.EQ.3)SR  =  <i. 03^(71 
IF(0RDER.EQ.'i)SR  =  2. 11596 
IF(0RDER.E0.5)SR=1 .50711 
IF(0RDER.E(J.6)SR  =  1. 25325 

ENDIF 

IF(NATTN.E().70)THEN 

IF(ORDER.EQ.2)SR=28.15950 
IF(0RDER.EQ.3)SR=5.87  251 
IF(0RDER.EQ.'i)SR  =  2.7<i7^9 
IF(0RDER.EQ.5)SR=1. 80680 
IF(0RDER.EQ.6)SR  =  l.<il799 


NEW05770 
NEW05780 
NEW05790 
NEM05800 
NEN05810 
NEW05820 
NEH05830 
NEW058'40 
NEH05850 
NEH0586  0 
NEW05870 
NEW05880 
NEW05890 
NEW05900 
NEW05910 
NEW05920 
NEW05930 
NEW059^0 
NEW05950 
NEW05960 
NEW05970 
NEW05980 
NEW05990 
NEW06000 
NEH06010 
NEW06020 
NEN06030 
NEW060^0 
NEW06050 
NEN06060 
NEH06  07  0 
NEW06080 
NEH06090 
NEH06100 
NEH06110 
NEW06120 
NEW06130 
NEH061<i0 
NEII06150 
NEH06160 
NEW06170 
NEH06180 
NEW06190 
NEW06200 
NEH06210 
NEW06220 
NEW06230 
NEW062^0 
NEW06250 
NEW06260 
NEII06270 
NEH06280 
NEH06290 
NEW06300 
HEH06310 
NEH06320 
NEH06330 
NEN06  3<iO 
NEH06350 
NE1I06360 
NEW06370 
NEH06  38  0 
NEH06390 
NEN06<400 
NEI'106410 
NEN06<<20 
NEII06'i30 
NEH06<*<«0 
NEW06<450 
NEW06460 
NEW06'i70 
NEH06<i80 
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FORTRAN 


17 


19 


ENDIF 
EMDIF 

SPI=3.1'H59 
RK=1.0/SR 
N=ORDER 
EN  =  N 

KPR=SQRT(1.0-RK»»2) 

QO=0.5K(1.0-SORT(KPR))/(1.0+S(3RT(KPR)) 
Q=QO+2.0XQOxx5+15.0*QOx»9+150.0XQO»xl3 
D=(10.0XX(0.1XASUBA)-1 . 0)/< 10 . OXK( 0 . 1«ASUBP)-1 . 0) 
LAM=(1.0/(2.0*EN))XALOG((10.0X«(0.05KASUBP)+1.0)/(10.0«»(0.05 

XASUBP)-1.0)) 
ISTOP=0 
DENSUM=0.0 
EHSUM=SINH(LAM) 
M=l 
EM=M 

ENUM=((-1.0)»XM)K(0K»(MX(M+1)))»SINH((2.0»EM+1.0)XLAM) 
DEH=2.0X((-1.0)»XM)«(QK*CM«X2))KCOSH(2.0XEMKLAM) 
ENSUM=ENSUM+ENUM 
DENSUM=DENSUM+DEN 
IF(M.GE.2)THEN 

RN=ABS(ENUM/ENSUM) 

RD=ABS(DEN/(1.0+DENSUM)) 

IF((RN.LE.1.0E-8).AND.(RD.LE.1.0E-8))THEN 
IST0P=1 

ENDIF 
ENDIF 
M=M+1 
IF(ISTOP.EQ.0)THEN 

GO  TO  17 
ENDIF 

SIGMAO=ABS((2.0X(QXX0.25)»ENSUM)/(1.0+DENSUM)) 
W=SQRT((1.0+RKX(SIGMAOXX2))X(1.0+(SIGMAOXX2/RK))) 
1N=M0D(N,2) 
IF(IN.EQ.O)THEN 

lR=N/2 
ELSE 

IR=(N-l)/2 
ENDIF 
DO  18  1=1, IR 

IF(ir).EQ.O)THEN 
EMU=I-0.5 

ELSE 
EMU  =  I 

ENDIF 

ISTOP=0 

DENSUM=0.0 

ENSUM=SIN(PIXEMU/EN) 

M  =  I 

EM  =  M 

ENUM=((-1.0)XXM)X(QXK(MXCM+1)))XS1N((2.0XEM+1 . 0)xPlxEMU/EN) 

DEN  =  2.0X((-1.0)XXM)X(Q?'.X(MXX2))XC0S(2.0KEM«PIKEMU/EN) 

ENSUM=E(ISUM+ENUM 

DENSUM=DENSUM+DEN 

IF(M.GE.2)THEN 

RN=ABS(ENUM/ENSUM)  ■' 

RD=ABS(DEN/(1.0+DENSUM)) 

IF((RN.LE.1.0E-8).AND.(RD.LE.1.0E-8))THEN 

ISTQP=1 
ENDIF 

ENDIF 

M  =  M+I 

IF(ISTOP.EQ.0)THEN 
GO  TO  19 

ENDIF 

OMEGA(I)=(2.0»(QXX0.25)*ENSUM)/(1. 

V(I)=SQRT((1.0-RKXOMEGA(I)XX2)X(1, 

A0(I)=1 .0/(OMEGA(I)XX2) 

DN  =  1  .  0  +  (SIGMAOXOMEGA(I))»(J(2 

B0(I)=((SIGMA0XV(I))»X2+{0MEGA(1)XW)KK2)/(DNKX2) 

B1(I)=(2.0XSIGMAOXV(I))/DN 


•O+DENSUM) 
.0-(0MEGA(I)xx2/RK))) 


NEN06<i90 
NEH06500 
NEH06510 
NEH06520 
NEW06530 
NEM065<«0 
NEW06550 
NEW06560 
NEW06570 
NEH06580 
NEW06590 
NEH06600 
NEM06610 
NEW06620 
NEM06630 
NEW066<tO 
NEW06650 
NEH06660 
NEW06670 
NEN06680 
NEW06690 
NEN06700 
NEW06710 
NEW06720 
NEW06730 
NEW067'iO 
NEH06750 
NEW06760 
NEW06770 
NEW06780 
NEW06790 
NEkl06800 
NEH06810 
NEI-106820 
NEW06830 
NEH068'iO 
NEH06850 
NEH06860 
NEW06870 
NEH06880 
NEN06890 
NEW06900 
NEN06910 
NEW06920 
NEH06930 
NEI>I069<<0 
NEW06950 
NEH06960 
NEH06970 
NEH06980 
NEN06990 
NEW07000 
NEW07010 
NEW07020 
NEW07030 
NEH070'iO 
NEH07050 
NEW07  06  0 
NEW07  07  0 
NEN07080 
NEH07090 
NEH07100 
HEW07110 
NEH07120 
NEW07130 
NEH071'iO 
NEW07150 
NEH07160 
NEW07170 
NEW07180 
NEW07190 
NEW07200 
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18      CONTINUE 
H0=1 .0 

IF(IN.EQ.O)THEN 
DO  20  1  =  1, IR 

H0=H0XB0CI)/A0(I) 

20  CONTINUE 
HO=HOX(10.0XK(-0.05»A5UBP)) 

ELSE 

DO  21  1=1, IR 

H0=H0XB0(I)/A0(I) 

21  CONTINUE 
H0=H0XSIGMA0 

ENDIF 

R=DBLE(SR) 
H0=DBLE(HO) 
A01  =  DBLE(AO(1).) 
B01=DBLE(BO(1)) 
B11=DBLE(B1(1)) 

I F( ( ORDER . EQ . 2 ) . OR . ( ORDER . EQ . 3) )THEN 
APROT(1,0)=HO*A01 
APROT(1,1)=O.OD+00 
APROT(l,2)=H0 
APROT(1,3)=O.OD+00 
IF(0RDER.E9.2)THEN 

APROT(2,0)=B01 

APR0T(2,1)=B11 

APR0T(2,2)=1.0D+00 
ENDIF 
IF(0RDER.EQ.3)THEN 

S0=DBLE(SIGMAO) 

APROT(2,0)=S0»B01 

APROTC2,1)=B01+(B11XSO) 

APROT(2,2)=SO+B11 

APROT{2,3)=1.0D+00 
ENDIF 
ENDIF 

IF((ORDER.EQ.«i)  .  OR  .  (ORDER  .  EQ  .  5)  )THEN 
A02=DBLE(AO(2)) 
B02=DBLE(Ea(2)) 
B12=DBLE(B1(2)) 
APROT(1,0)=HOXA01)(A02 
APROT(1,1)=O.OD+00 
APROT(1,2)=HO*(A01+A02) 
APROT(l,3)=0.0D+00 
APROTCl,'i)=H0 
APROT(l,5)=0.0D+O0 
IF(ORDER.EQ.'i)THEN 

APROT(2,0)=B01XB02 

APROT(2,1)  =  (B11)«B02)  +  (B01»B12) 

APROT(2,2)=B01+B02+(B11XB12) 

APR0T(2,3)=B11+B12 

APR0T(2,'4)  =  1.0D+00 
ENDIF 
IF(0RDER.EQ.5)THEN 

S0=DBLE(SIGMAO) 

APROT(2,0)=SOXB01KB02 

APROTC 2, l)=(S05<((BllxB02)+( 601*612 )))+(B01XB02) 

APROT(2,2)  =  (SO)((601  +  B02+(B1IX612)))  +  (611)(B02)  +  (B01X612) 

APROT(2,3)=(S0)((Bll+B12))+B01+B02+CeilxB12) 

APR0T(2,<4)=SO  +  B11  +  B12 

APROT(2,5)=1.0D+00 
ENDIF 
ENDIF 

1F((0RDER.EQ.6).0R.(0RDER.EQ.7))THEN 
A02=DBLE(AO(2)) 
A03  =  DBLE(AO(3))  i;» 

B02=DBLE(BO(2)) 
B03=DBLE(BO(3)) 
612=DBLE(61(2)) 
B13=DBLE(B1(3)) 
APROT(1,0)=HOKA01XA02XA03 
APROT(1,1)=0,OD+00 


NEH07210 
NEW07220 
NEW07230 
NEW072<iO 
NEW07250 
NEW07260 
NEW07270 
NEH07280 
NEW07290 
NEW07300 
NEH07310 
NEH07320 
NEW07330 
NEW073<40 
NEW07350 
NEW07360 
NEW07370 
NEW07380 
NEW07390 
NEW07<»00 
NEW07<il0 
NEW07<i20 
NEW07<i30 
NEH07^<iO 
NEW07<»50 
NEH07^60 
NEH07<i70 
NEK07480 
NEW07^90 
NEH07500 
NEH07510 
NEH07520 
NEW07530 
NEH075^0 
NEH07550 
NEH07560 
NEH07570 
NEH07580 
NEW07590 
NEH07600 
NEW07610 
NEH07620 
NEH07630 
NEll076<iO 
NEl')07650 
NEN07660 
NEW07670 
NEH07680 
NEII07690 
NEW07700 
NEH07710 
NEH07720 
NEH07730 
NEWa77<i0 
NEW07750 
NEW07760 
NEW07770 
NEH07780 
NEN07790 
NEII07800 
NEW07810 
NEH07820 
NEW07830 
NEH078<(0 
NEH07850 
NEH07860 
NEW07870 
NEH07880 
NEH07890 
NEN07900 
NEN07910 
NEW07920 
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26 
25 

C 
C 


50 


51 

C 
C 
C 
C 
C 

c 


61 

60 


APR0T(1,2)=H0»((A02KA03)+(A01K(A02+A03))) 
APROT(l,3)=0.0D+00 
APROT(1,^)=HOX(A01+A02+A03) 
APROT(l,5)=0.0D+00 
APROT(l,6)=H0 

I FC  ( ORDER  .  EQ  .  6  )  .  OR  .  ( ORDER .  EQ .  7  )  )THEN  -- 
APROT(2,0)=B01KB02»B03 

APROT(2,1)=(B11*B02«B03)+(B01»((B02»B13)+(B12)«B03))) 
APROT(2,2)=(B02*B03)+(B11«((B02KB13)+{B12XB03)))+ 
:        (B01X((B02+B03)+(B12XB13))) 

APROT(2,3)=B01X(B12+B13)+(B11«((B02+B03)+(B12XB13)))+ 
C        ((B02KB13)+(B12»B03)) 

APROT(2,<i)=B01+(Bll«(B12+B13))+((B02+B03)+(B12«B13)) 
APR0T(2,5)=B11+B12+B13 
APROT(2,6)=1.0D+00 
ENDIF 

IF(aRDER.EQ.7)THEN 
SO=DBLE(SIGMAO) 
APROT(2,7)=1.0D+00 

APROT(2,6)=(S0XAPROT(2,6))+APROT(2,5) 
APR0T(2,5)=(S0XAPROT(2,5))+APR0T(2,<4) 
APROTC2,<i)  =  (S0xAPROT(2,^))+APROT(2,3) 
APR0T(2,3)=(S0XAPR0T(2,3))+APR0T(2,2) 
APROT(2,2)=(S0XAPROT(2,2))+APROT(2,l) 
APR0T(2,l)=(S0xAPROT(2,l))+APROT(2,0) 
APROT(2,0)=SOXAPROT(2,0) 
ENDIF 
ENDIF 
ENDIF 

NORMALIZE  THE  ELLIPTIC  PROTOTYPES 
IF(TYPE.EQ.5)THEN 
DO  25  1=1,2 
DO  26  J=0, ORDER 

APROT(I,J)=APROT(I,J)/((DSQRT(R))KXJ) 
CONTINUE 
CONTINUE 
ENDIF 
APR0T(TYPE=6)  ALLOHS  THE  INPUT  OF  ANY  TYPE  OF  PROTOTYPE  Fill 

BY  INPUTING  THE  COEFFICIENTS  OF  S 
IF(TYPE.EQ.6)THEN 
WRITE(6,K) 

•THIS  OPTION  ALLOWS  YOU  TO  USE  ANY  ANALOG   • 
•  PROTOTYPE  OF  ORDER  1-10        • 
•»K  COEFFICIENTS  MUST  BE  IN  DOUBLE  PRECISION  • 
'     FORMAT  (0.123<456789D+00)  XX  • 


'INPUT  THE  NUMERATOR  COEFFICIENTS  FIRST' 


CCLRSCRN') 
'NOW  INPUT  THE 


DENOMINATOR  COEFFICIENTS' 


HRITE(6,x) 
WRITE(6,x) 
WRITE(6,x) 
WRITE(6,K) 
WRITE(6,X) 
WRITE(6,x) 
DO  50  1=0, ORDER 

WRITE(6,x)  'WHAT  IS  THE  COEFFICIENT  OF  SXX',I,'?' 

READC5,X)  APR0T(1,I) 
CONTINUE 
WRITE(6,x) 
CALL  EXCMS 
WRITE(6,») 
WRITE(6,») 
DO  51  J=0, ORDER 

WRITE(6,x)  'WHAT  IS  THE  COEFFICIENT  OF  SXX«,J,'T' 

READ(5,X)  APR0T(2,J) 
CONTINUE 
ENDIF 

XXXXXSOfXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXX 
XXXXXXXKX^XXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DEFINE  ELEMENTS  OF  THE  BILINEAR  TRANSFORMATION  MATRIX 

XXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

FOR  ALL  THE  BILINEAR  TRANSFORATION  MATRICES,  THE  FIRST  COLUMN 
HAS  ALTERNATING  +1'S  AND  -I'S.  THE  LAST  COLUMN  HAS  ALL  +1'S 
DO  60  1=1, ORDER 
DO  61  J=0,I 
BILINR(I,J,0)=(-1.0)XXJ 
BILINR(I,J,I)=1 .0 
CONTINUE 
CONTINUE 


NEW07930 
NEW079<40 
NEW07950 
NEW07960 
NEW07970 
NEH07980 
NEW07990 
NEM08000 
NEN08010 
NEW08020 
NEW08030 
NEW080<40 
NEH08050 
NEW08060 
NEW08070 
NEW08080 
NEW08090 
NEW08100 
NEW08110 
NEW08120 
NEH08130 
NEW081<)0 
NEW08150 
NEW08160 
NEW08170 
NEW08180 
NEW08190 
NEW08200 
NEW08210 
NEW08220 
NEW08230 
NEW082^0 
NEW08250 
NEW08260 
NEW08270 
NEW08280 
NEW08290 
NEW08500 
NEW08310 
NEW08320 
NEW08330 
NEW083<iO 
NEW08350 
NEW08360 
NEW08370 
NEW08380 
NEH08390 
NEW08400 
NEW08<4lO 
NEW08<i20 
NEW08'i30 
NEH08'(^0 
NEW08'i50 
NEW08<460 
NEW08<i70 
NEW08<«80 
NE1I08^90 
NEW08500 
NEW08510 
NEW08520 
NEW08530 
NEW085<*0 
NEW08550 
NEW08560 
NEW08570 
NEW08580 
NEW08590 
NEW08600 
NEW08610 
NEW08620 
NEW08630 
NEH086<«0 


?» 
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C 
C 

c 


77 
76 
C 


78 
75 
C 
C 


81 
80 
C 
C 

c 


102 

101 
100 
C 

c 


106 
105 


111 


112 


C 
C 
C 

c 
c 
c 


NEW08650 

FIRST  ORDER  MATRIX  HAS  ALREADY  BEEN  LOADED  CONTINUE  INTERATION    NEW08660 

UNTIL  YOU  HAVE  THE  APPROPRIATE  MATRIX  NEW08670 

DO  75  K=2, ORDER  NEW08680 

KM1=K-1  NEW08690 

LOAD  ELEMENTS  IN  INTERIOR  COLUMNS,  ALL  ROWS  EXCEPT  LAST  ONE  NEW08700 

DO  76  L=0,KM1  NEW08710 

DO  77  M=1,KM1  NEW08720 

MM1=M-1  NEW08730 

BILINR(K,L,M)=BILINR(KM1,L,M)+BILINR(KM1,L,MM1)  NEW087^0 

CONTINUE  NEW08750 

CONTINUE  NEW08760 

NOW  LOAD  THE  LAST  ROW  NEW08770 

DO  78  N=1,KM1  NEW08780 

BILINR(K,K,N)=BILINR(K,0,N)X(X-1.0D+00)««(K+N))  NEW087  90 

CONTINUE  NEW08800 

CONTINUE  NEW08810 

XXXXXXXXKXKXXXXXXiCXXXXXXXXXX^CXKXKKXXXXiOCXXXKXXXXKMXKXKKXXXXXXXKX   NEW08820 

THE  FOLLOWING  CLEARS  DPROT  PRIOR  TO  LOADING  NEW08830 

DO  80  K  =  l,2  NEW088'i0 

DO  81  L=0,10  NEW08850 

DPROT(K,L)=0.0  NEW08860 

CONTiriUE  NEW08870 

CONTINUE  NEW08880 
XXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXX  NEW08890 
DIGITIZE  THE  ANALOG  PROTOTYPE  USING  THE  BILINEAR  TRANSFORMATION   NEW08900 

XXXXXXXXXXXXXXXXXXXXXXJ(XXXXXXXXXXXXXXXXXXXXXXXXXKXXXKXXXXXKXXXXX   NEW08910 

DO  100  N=l,2  NEN08920 

DO  101  J=0, ORDER  NEW08930 

SUMR  =  0.0  NEH089't0 

DO  102  1=0, ORDER  NEH08950 

TEMP=APROT(N,I)xBILINR(ORDER,I,J)  NEW0896  0 

SUMR=SUMR+TEMP  NEW08970 

CONTINUE  NEH08980 

DPROTCN, J)=SUMR  NEW08990 

CONTINUE  NEW09000 

CONTINUE  NEW09010 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  NEH09020 

OPTION  TO  SHOW  DIGITAL  PROTOTYPE  COEFFICIENTS  NEH09030 

CALL  EXCMS  CCLRSCRN')  NEW090<i0 

HRITE(6,x)  'THE  DIGITAL  PROTOTYPE  HAS  BEEN  COMPUTED.  DO  YOU'  NEII09050 

HRITE(6,x)  'WISH  TO  CHECK  THE  PROTOTYPE  COEFFICIENTS? •  NEW09060 

WRITE(6,x)  •   XX  1-YES/2-N0  xx'  NEH09070 

READ(5,x)  LOOK  NEW09080 

IFCLOOK.EQ.DTHEN  NEH09090 

MAKE  DENOMINATOR  COEFFICIENT  OF  ZXXORDER  =1.0  NEW09100 

DO  105  K=l,2  NEW09110 

DO  106  L=0, ORDER  NEH09120 

DPROT ( K, L)= DPROT (K,L)/DPROT( 2, ORDER)  NEW09130 

CONTINUE  NEH091<iO 

CONTINUE  NEW09150 

CALL  EXCMS  CCLRSCRN')  NEH09160 

WRITE(6,»)  'THE  DIGITAL  PROTOTYPE  NUMERATOR  COEFFICIENTS  ARE>'    NEH09170 

DO  111  N=ORDER,0,-1  NEW09180 

WRITE(6,902)  DPR0T(1,N),N  NEH09190 

CONTINUE  NEW09200 

HRITE(6,K)  NEH09210 

HRITE(6,x)  'THE  DIGITAL  PROTOTYPE  DENOMINATOR  COEFFICIENTS  ARE.'  NEH09220 

DO  112  M=ORDER,0,-1  NEH09230 

WRITE(6,902)  DPR0T(2,M),M  NEW092^0 

CONTINUE  NEW09250 

IIRITE(6,x)  NEH09260 

WRITE(6,X)  'XX  1  TO  CONTINUE  XX'  NEW09270 

READ(5,x)  LOOK  NEW09280 

ENDIF  NEW09290 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXX   NEW09300 

NOW  DO  THE  DIGITAL  TO  DIGITAL  TRANSFORMATIONS  NEW09310 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXK   NEW09320 

IF((PASS.EQ.l) .0R.(PASS.EQ.2))THEN  NEW09330 

DEFINE  ELEMENTS  OF  THE  LOWPASS  AND  HIGHPASS  DIGITAL-DIGITAL  NEW093<i0 

TRANSFORMATION  MATRIX  NEW09350 

DEFINE  THE  ELEMENTS  FOR  FIRST  AND  LAST  COLUMNS  NEW09360 
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171 

170 

C 

C 


177 
176 


178 

175 
C 

c 
c 


181 
180 


187 

186 
185 

C 
C 

C 
C 


226 
225 
C 


DO  170  1=1,10  NEW09370 

DO  171  J=0,I  NEW09380 

LPHPTR(I,J,0)  =  (((-1.0)>«ALPHA)i<XJ)  NEW09390 

IMJ  =  I-J  .,  NEW09<i00 

LPHPTR(I,IMJ,I)  =  LPHPTRCI,J,0)  NEW09<ilO 

CONTINUE  NEW09<»20 

CONTINUE  NEW09<i30 

FIRST  ORDER  MATRIX  HAS  HF^.H    DEFINED  ABOVE  NEW09^<iO 

SECOND  TO  TENTH  ORDER  BELOW  NEW09<)50 

DO  175  K  =  2,10  NEW09'i60 

KM1  =  K-1  NEW09<470 

DO  176  L  =  0,KM1  NEW09<i80 

DO  177  M=1,KM1  NEW09<)90 

MM1=M-1  NEW09500 

LPHPTR(K,L,M)=LPHPTR(KMl,L,M)+((LPHPTR(KMl,L,MMl))K(-1.0D+00)  NEW09510 

CKALPHA)  NEW09520 

CONTINUE  NEW09530 

CONTINUE  NEW095^0 

DO  178  N=1,KM1  NEW09550 

KMN=K-N  NEW09560 

LPHPTR(K,K,N)=LPHPTR(K,0,KMN)  NEH0957  0 

CONTINUE  NEW09580 

CONTINUE  NEH09590 

XXXXKXXXXXXK^XKitXXXKiOiXKKXXXX^CXXXXXXXKKXXXXKXXXXKKKXKXXXXXKXKXXX   NEH096  00 

CALCULATE  THE  FINAL  FILTER  NEW09610 

FOR  HIGHPASS,  THE  PROTOTYPE  MUST  SUBSTITUTE  (-Z)  FOR  (Z)  NEH09620 

IF(PASS.EQ.2)THEN  NEH09630 

DO  180  K  =  l,2  NEW096<«0 

DO  181  L=0, ORDER  NEH09650 

DPROT(K,L)=DPROT(K,L)X((-1.0)XXL)  NEW0966  0 

CONTINUE  NEW09670 

CONTINUE  NEH09680 

ENDIF  NEW09690 

DO  185  K  =  l,2  NEtl09700 

DO  186  L=0, ORDER  NEH097I0 

SUMR=0.0D+00  NEW09720 

DO  187  M=0, ORDER  NEW09730 

TEMP=DPROT{K,M)XLPHPTR( ORDER, M,L)  NEH097^0 

SUMR=5UMR+TEMP  NEW09750 

CONTINUE  NEH09760 

DFLTR(K,L)=SUMR  NEW09770 

CONTINUE  NEW09780 

CONTINUE  NEW09790 

ENDIF  NEH09800 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  NEW09810 

BANDPASS  AND  BANDSTOP  FORMULAS  NEH09820 

IF((PASS.E0.3)  .OR.(PASS.EQ.<i))THEN  NEH09830 

LOAD  THE  BANDPASS/BANDSTOP  TRANSFORMATION  MATRIX  NEH098<40 

FIRST  LOAD  THE  FIRST  AND  LAST  COLUMNS  NEW09850 

DO  225  1=1,10  NEH09860 

DO  226  J=0,I  NEH09870 

BPBSTRCI, J,0)=CXXJ  NEW09880 

12=2X1  NEH09890 

IMJ  =  I-J  NEII09900 

BPBSTR(I,IMJ,I2)=BPBSTR(I,J,0)  NE1109910 

CONTINUE  NEH09920 

CONTINUE  NEH09930 

LOAD  THE  REST  OF  THE  FIRST  ORDER  MATRIX  NEW099^0 

BPBSTR(1,0,1)=(-1 .0D+00)XB  NEW09950 

BPBSTR(1,1,1)=(-1.0D+00)XB  NEH0996  0 

LOAD  THE  2ND-10TH  ORDER  TRANSFORMATION  MATRICES  NEH09970 

DO  230  K=2,I0  NEW09980 

K2  =  2»'K  NEW09990 

K2M1=K2-1  NEHIOOOO 

KM1=K-1  NEWIOOIO 

K2M2=K2-2  NEW10020 

KM12=2XKM1  NEW10030 

KM12Mi=KM12-l  NEW100<»0 

LOAD  THE  SECOND  AND  SECOND  TO  LAST  COLUMNS  NEW10050 

DO  231  L  =  0,KM1  NEm0060 

BPBSTR(K,L,1)=BPBSTR(KM1,L,1)+(BPBSTR(KM1,L,0)X(-1.0D+00)»B)    NEW1007  0 

BPBSTR(K,0,K2M1)=(BPBSTR(KM1,0,KM12)X(-1.0D+00)XB)+  NEW10080 
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231 

C 


237 
236 
C 


238 
230 
C 
C 

C 


211 
210 


217 

216 
215 

C 
C 

c 


251 
250 


300 


301 


C      (BPBSTR(KM1,0,KM12M1)»C) 
KML=K-L 

BPBSTR(K,KML,K2M1)=BPBSTR(K,L,1) 
CONTINUE 
NOW  LOAD  INTERIOR  ELEMENTS  EXCEPT  FOR  LAST  ROW 
DO  236  L=0,KM1 
DO  237  M=2,K2M2 
MM1=M-1 
MM2=M-2 

BPBSTR(K,L,M)=BPBSTR(KM1,L,M)+(BPBSTR(KM1,L,MM1)X(-1.0D+00)« 
C      B)+(BPBSTR(KM1,L,MM2)KC) 
CONTINUE 
CONTINUE 
NOW  FILL  THE  LAST  ROW 
DO  238  N=1,K2MI 
K2MN=K2-N 

BPBSTR(K,K,K2MN)=BPBSTR(K,0,N) 
CONTINUE 
CONTINUE 

XXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXKKXXXKXXiCKXXKXXKXKXXXXXXXXX 

CALCULATE  THE  FINAL  FIITFR"' 

FOR  THE  BANDPASS  CASE  REPLACE  (Z)  WITH  (-2) 
IF(PASS.EQ.3)THEN 
DO  210  M=l,2 
DO  211  N=0, ORDER 

DPROT(M,N)=DPROT(M,N)X((-1.0D+00)XXN) 
CONTINUE 
CONTINUE 
ENDIF 
0RDER2=2X0RDER 
DO  215  K=l,2 
DO  216  L=0,ORDER2 
SUMR-O.OD+00 
DO  217  M=0, ORDER 
TEMP=DPROT(K,M)XBPBSTR( ORDER, M,L) 
SUMR=SUMR+TEMP 
CONTINUE 
DFLTR(K,L)=SUMR 
CONTINUE 
CONTINUE 
ENDIF 

XXXXXXXXSfXXXXXXX 

OUTPUT  OPTIONS 

KXXXXXXXXXXXXXXX 

IF((PASS.EQ.3).0R.(PASS.EQ.1))THEN 

END=2X0RDER 
ENDIF 
IF((PASS.EQ.1).0R.(PASS.EQ.2))THEN 

END=ORDER 
ENDIF 

MAKE  DENOMINATOR  COEFFICIENT  OF  ZXXEND  EQUAL  TO  1 
DO  250  M=l,2 

DO  251  N=0,END 
DFLTR(M,N)=DFLTR(M,N)/DFLTR(2,END) 

CONTINUE 
CONTINUE 

CALL    EXCMS    CCLRSCRN') 
HRITE(6,x) 
HRITE(6,xj 
WRITE(6,x) 

HRITE(6,x)    •  NUMERATOR   COEFFICIENTS. • 

DO    500    I=END,0,-1 

HRITE(6,902)    DFLTR(1,I),1 
CONTINUE 

WRITE(6,x)  «     DENOMINATOR  COEFFICIENTS! • 
DO  301  J=END,0,-1 

WRITE(6,902)  DFLTR(2,J),J 
CONTINUE 
WRITE(6,x) 

WRITE(6,x)  •  XX  ENTER  1  TO  CONTINUE  XX' 
READ(5,x)  LOOK 
XXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXKXXXXXXXXXXXX 


•THE  FINAL  DIGITAL  FILTER  IS:' 


NEW10090 
NEWIOIOO 
NEWlOllO 
NEW10120 
NEH10130 
NEWlOllO 
NEW10150 
NEW10160 
NEW10170 
NEW10I80 
NEW10190 
NEW10200 
NEW10210 
NEW10220 
NEW10230 
NEW10210 
NEW10250 
NEW10260 
NEW10270 
NEW10280 
NEW10290 
NEH10300 
NEH10310 
NEW10320 
NEW10330 
NEH10310 
NEH10350 
NEW10360 
NEW10370 
NEII10380 
NEW10390 
NEWIOIOO 
NEWlOllO 
NEW10120 
NEW10130 
NEWlOllO 
NEW10150 
NE111016  0 
NEW10170 
NEW10180 
NEW10190 
NEW10500 
NEW10510 
NEW10520 
NEH10530 
NEH10510 
NEW10550 
NEW10560 
NEW10570 
NEW1058  0 
NEW10590 
NEW10600 
NE1I10610 
NEIH0620 
NEW10630 
NEW10610 
NEN10650 
NEH10660 
NEH106  7  0 
NEII10680 
NEH106  9  0 
NEW10700 
NEW10710 
NEW10720 
NEW10730 
NEW10710 
NEW10750 
NEW10760 
NEW10770 
NEW10780 
NEW10790 
NEW10800 
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310 


320 


OPTION  TO  SPOT  CHECK  FILTER  OUTPUT  MAGNITUDE  AND  FREQUENCY 
CALL  EXCMS  CCLRSCRN') 

HRITE(6,»)  'THE  TRANSFER  FUNCTION  FOR  THIS  FILTER  HAS  BEEN' 
WRITE(6,K)  'COMPUTED.  DO  YOU  WANT  TO  SPOT  CHECK  THE  MAGNITUDE' 
WRITE(6,x)  'AND  PHASE  OF  THE  FILTER  OUTPUT  FOR  A  GIVEN  ' 
HRITE(6,K)  'FREQUENCY?' 
WRITE(6,«)  '     K«  1-YES/2-N0  XX  ' 
READ(5,x)   LOOK 
IFCLOOK.EQ.DTHEN 
CALL  EXCMS  CCLRSCRN') 
IF(AF0RDF.EQ.1)THEN 

WRITE(6,X)  'WHAT  ANALOG  FREQUENCYCIN  KHZ)' 
WRITE(6,x)  •  DO  YOU  WANT  TO  CHECKT ' 
WRITE(6,J()  •    XK  DECIMAL  XX" 
READ(5,X)   F 
DF=(F/SSFREQ)X2.0XSPI 
ENDIF 
IF(AF0RDF.E0.2)THEN 

WRITE(6,x)  'WHAT  DIGITAL  FREQUENCY  DO  YOU  WANT  TO  CHECKT' 
WRITE(6,x)  •   XX  DECIMAL  XX  ' 
READ(5,X)  DF 
ENDIF 
X=COS(DF) 
Y=SINCDF) 
Z=CMPLX(X,Y) 
XN=REAL(DFLTR(1,0)) 
YN=0.0 

XD=REAL(DFLTR(2,0)) 
YD=0.0 

HOFZH=CMPLX(XN,YN) 
HOFZD=CMPLX(XD,YD) 
DO    320    N  =  l,ErtD 
ZN=REAL(DFLTR(l,N))x(Z«*N) 
ZD=REAL(DFLTR(2,N))K(Z««N) 
HOFZN=HOFZN+ZN 
HOFZD=HOFZD+ZD 
COHTINUE 

HOFZ=:HOFZN/HOFZD 

FMAG=SQRT( ( REAL ( H0FZ)KX2 )+( AIMAG( H0FZ)«X2) ) 
FMAGDB=20.0*ALOG10(FMAG) 
FPHSE=ATAN(AIMAG(HOFZ)/REAL(HOFZ)) 
WRITE(6,x)  'THE  FILTER  OUTPUT  WILL  BE:' 
WRITE(6,91<i)  FMAGDB 
WRITE(6,915)  FPHSE 
WRITE(6,x) 

WRITE(6,*)    'DO   YOU   WANT    TO   CHECK    ANOTHER    FREQUENCY?' 
WRITE(6,K)    '  XX    1-YES/2-N0    xx    • 

READ(5,x)    LOOK 
IFCLOOK.EQ.DGO   TO    310 
ENDIF 

XitXKXJtXXXXXXXX^tXXXKKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXX 

OPTION  TO  LOAD  FREQUENCY  RESPONSE  INTO  A  DATA  FILE 

CALL  EXCMS  CCLRSCRN') 

'DO  YOU  WANT  THE  FILTERS  FREQUENCY  RESPONSE  ' 
'DATA  TO  BE  ENTERED  INTO  A  DATA  FILE  FOR  USE' 
'BY  A  PLOTTING  PROGRAM?  • 

'  NOTE:  THE  FILE  NAME  WILL  BE  FILTER,  FILE  TYPE' 
'      WILL  BE  DATA,  100  POINTS  WILL  BE  CALCULATED' 

•  FOR  A  DIGITAL  FREQUENCY  FROM  0  TO  PI,  THE  FIRST' 
'      COLUMN  WILL  BE  THE  DIGITAL  FREQUENCY  (EXPRESSED 
'      AS  A  FRACTION  OF  THE  SAMPLING  FREQUENCY  0-1/2)  ' 

•  THE  SECOND  COLUMN  WILL  BE  THE  MAGNITUDE,  AND' 
'      THE  THIRD  COLUMN  WILL  BE  THE  PHASEC DEGREES) ' 
'      XX  1-MAKE  DATA  FILE/2-D0N"T  BOTHER  xx' 


WRITE(6,x) 

WRITE(6,x) 

WRITE(6,x) 

HRITE(6,x) 

WRITE(6,») 

WRITE(6,*) 

WRITE(6,x) 

WRITE(6,X) 

WRITE(6,x) 

WRITE(6,K) 

WRITE(6,X) 

READ(5,x) 

IFCLOOK.EQ 


LOOK 
DTHEN 
DO  350  1=0,100 
DF=REAL(I)X(5. 0/1000.0) 
DFR=REAL(I)X(SPI/100.0) 
X=COS(DFR) 
Y=SIN(DFR) 
Z=CMPLX(X,Y) 


NEW10810 
NEW10820 
NEW10830 
NEW108<iO 
NEW10850 
NEW10860 
NEW10870 
NEW10880 
NEW10890 
NEW10900 
NEW10910 
NEW10920 
NEW10930 
NEW109<i0 
NEW10950 
NEW10960 
NEW10970 
NEW10980 
NEW10990 
NEWUOOO 
NEWnOlO 
NEW11020 
NEW11030 
NEWllO^iO 
NEW11050 
NEW1106  0 
NEW11070 
NEW11080 
NEW11090 
NEWlllOO 
NEWllllO 
NEW11120 
NEW11130 
NEWlllOO 
NEW11150 
NEW11160 
NEW11170 
NEW11180 
NEW11190 
NEW11200 
NEW11210 
NEW11220 
NEW11230 
NEW112aO 
NEW11250 
NEW1126  0 
NEW11270 
NEtlll280 
NEW11290 
NEW11300 
NEW11310 
NEW11320 
NEW11330 
NEW113^0 
NEW11350 
NEW1136  0 
NEW1137  0 
NEtlll38  0 
NEW11390 

'NEWlllOO 
NEWll'ilO 
NEWll<i20 
NEW11^30 

. NEWll^^O 
NEW11<J50 
NEWll<i60 
NEWll<i70 
NEWll<t80 
NEW11^90 
NEW11500 
NEW11510 
NEH11520 
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XN=REAL(DFLTR(1,0))  NEN11530 

YN  =  0.0  NEWllS'iO 

XD=REAL(DFLTR(2,0))  NEW11550 

YD=0.0  NEH1156  0 

HOFZN=CMPLX(XN,YN)  NEW11570 

HOFZD=CMPLX(XD,YD)  NEW11580 

DO  351  N=1,END  NEW11590 

ZN=REAL(DFLTR(1,N))KCZXKN)  NEW116  00 

ZD=REAL{DFLTR(2,N))K(ZXXN)  NEW11610 

HOFZN=HOFZN+ZN  NEW11620 

HOFZD=HOFZD+ZD  NEW11630 

351      CONTINUE  NEW116<«0 

HOFZ^HOFZN/HCFZD  NEW11650 

FMAG=S(3RT((REAL(HOFZ)XK2)  +  (AIMAG(HOFZ)XX2))  NEmi66  0 

IF(REAL(HOFZ).NE.0.0)THEN  NEW1167  0 

FPHSER=ATAN(AIMAG(HOFZ)/REAL(HOFZ))  NEH11680 

ELSEIFCREALCHOFZ)  .EQ.O.OTHEN  NEW11690 

IF(AIMAG(HOFZ).GT.O.O)FPHSER=1. 570796  NEmi700 

IF(AIMAG(HOFZ).LT.0.0)FPHSER=-1.57  0796  NEW11710 

IF(AIMAG(HOFZ).EQ.O.O)FPHSER=0.0  NEW11720 

ENDIF  NEW11730 

FPHSE  =  FPHSERS((18  0.0/SPI)  NEW117<iO 

NRITE(8,916)  DF, FMAG, FPHSE  NEW11750 

350     CONTINUE  NEW11760 

ENDIF  NEW11770 

C       KXXXKXXKXXXKKXXXXXXXXXKXKXXXXXXXXXKKXKXXKKXXXXKXXXXXXXXKXXXXXXXK   NEW11780 

901  FORMAT(5X,F10.5, •  Sxx',11)  NEW11790 

902  FORriAT(5X,F10.5,  '  Zxx',11)  NEW11800 

903  F0RMAT(3X,  •ALPHA=  \F9.<i)  NEWIISIO 

904  F0RMAT(3X, 'K^  •,F9.<i)  NEH11820 

905  F0RMAT(3X,'THE  LOHPASS  DIGITAL  PROTOTYPE  FREQ.  IS  ',F5.i)  NEW11830 

914  F0RMAT(3X, 'MAGNITUDE  =   ',F8.2,'  DB  •  )  NEH11840 

915  F0RMAT(3X, 'PHASE  =   ',F6.2,'  RADIANS')  NEW11850 

916  FORMAT(3X,F4.3,10X,F5.3,10X,F6.2)  NEH11860 
STOP  NEH11870 
END  NEWllSSO 
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